Gut microbiota maturation and early behavioral and cognitive development

Abstract

The presence of gut microbiota-brain-axis has been widely reported. However, few studies have focused on uncovering the potential associations during a time-period that our brain and gut microbiota undergo rapid maturation. We evaluated the potential associations between characteristics of gut microbiota and cognition and temperament using an accelerated longitudinal design in typically developing children over 0–3 years of age. Specifically, we extracted gut microbiota characteristics at three scale levels: diversity measures, microbial networks, and subject-wise longitudinal trajectory features, shedding light on how associations between cognition/temperament and gut microbiota may differ at global (diversity), ecological (microbial networks) and subject-wise levels. Our findings illustrated that associations between gut microbiota and temperament/cognition varied with the analytical approaches and highlighted differential gut microbial features in association with cognition and temperament traits—diversity measures and microbial networks largely with cognition while subject-wise trajectories with temperament. In addition, Ruminococcus bromii exhibited significant associations with cognitions spanning over multiple subdomains. Finally, the associations of gut microbiota with temperament and cognition converge on the potential interplay of language ability and social behaviors and highlight the importance of age-appropriate gut microbiota on early cognition/temperament development.

Introduction

Several lines of evidence from both human and animal studies have supported the presence of the gut-brain-axis—bidirectional interactions between the brain and gut microbiome^{1, 2, 3–4}. Diverse host-bacteria co-metabolites and hormones, like 5-HT, released by the enteroendocrine cells and effects of the vagus nerve could potentially affect brain health^3,5. In contrast, stress-related hormones released by the brain and mucosal immune components could alter gut microbial communities with reciprocal effects on brain functions^1,5. Causal relationships between gut microbiome maturation and diverse physiological measures and health conditions have been implicated, which highlight the importance of an age-appropriate gut microbiome maturation^1,6,7. Furthermore, numerous studies have reported potential relations between gut microbiota and symptoms of autistic children, particularly in the families of Bacteroidaceae,, Ruminococcaceae, and Lachnospiraceae^{8, 9, 10–11}. Finally, Johnson¹² found that subjects with larger social networks have a higher gut microbiota diversity, suggesting a possible link between gut microbiota and personality traits, which were further reinforced by a seminal review by Sherwin et al¹³ advocating reciprocal links between the gut microbiota and the social brain¹⁴. Collectively, current evidence has signified the potential interplay between the gut microbiome, the host’s general health, and social behavior and personality traits. However, most studies have focused on adults and young children with relatively few studies uncovering how gut microbiota maturation processes may be associated with brain development during infancy.

The gut microbiome undergoes a rapid increase of diversity, particularly between the ages of 4–12 months. At age 3, age specific community types begin to form, and maturation towards an adult like composition during the first five years of life¹⁵. Importantly, the temporal dynamic of gut microbiota maturation strikingly parallels rapid brain development during infancy^1,4,16. The first years of life represent one of the most critical periods of brain development^1,4. The brain undergoes rapid development with increasing brain volumes, myelination of white matter¹⁷, rapid maturation of basic and high order cognition¹⁸, and emergence of temperament. Temperament, referring to personality traits, has been alluded to as a pattern of emotion and cognition behaviors and assessed during childhood has been shown to predict personality and risk for social-emotional problems later in life¹⁹ such as anxiety²⁰, depression²¹ and ADHD²². Furthermore, Izard et al. reported significant relationships between cognitive ability and temperament in 7 years old children from economically disadvantage families²³, who arguably are more likely to experience dysbiosis due to environmental conditions and diet²⁴. Nevertheless, results in the literature have been conflicting; increased alpha diversity was associated with reduced fear reactivity and negative emotions²⁵ but also with lower cognitive ability and language functions²⁶, increased phylogenic diversity with higher surgency scores²⁷, and a decrease in Prevotella with increased behavioral problems²⁸. Discrepancies in the study designs, age ranges of study cohorts, geographic locations, and adjustments of different covariates may have accounted for the contradictory results. In addition, these aforementioned studies consisted of several limitations. Specifically, 16 s rRNA was employed by most studies, making it difficult to gain insights into the effects of microbes at the species level. Considering the close relation between cognition and temperament²⁹, further studies to simultaneously evaluate the potential associations between gut microbiota and cognition/temperament and how the associations may differ between cognition and temperament are needed. Few studies employed a longitudinal design, making it difficult to assess features of gut microbiota maturation trajectories at an individual level. It is highly plausible that features of gut microbiota developmental trajectories could yield important insights into development of cognition and temperament at the individual level. Finally, the gut microbiota is well recognized as an ecosystem where extensive interactions among species of microbes exist, underscoring the importance of considering features characterizing the interactions of microbes in the gut microbial community in relation to cognition and temperament.

With an accelerated longitudinal design, fecal samples were collected from typically developing children 0–3 years of age and analyzed using whole-genome shotgun methods, we aimed to uncover potential associations between cognition/temperament and characteristic of gut microbiota. We first considered how microbial compositions of gut microbiota at a global level may be associated with cognition and temperament using the widely used community diversity measures. While the diversity measures provide insights into the overall compositions of gut microbes, they do not capture potential interactions among microbes. It has been widely documented that microbes in the gut could either competitively or synergistically interact with each other³⁰, which in turn may affect their host’s physiology. In particular, the gut microbiota maturation has been characterized by sequential progression of different community structures during infancy^15,31. Such characteristics cannot be captured fully by microbiota diversity measures. To this end, network-based approaches taking into account the potential interactions among gut microbes were employed, enabling direct assessments of how characteristics of gut microbial networks may be associated with cognition and temperament. Finally, it has been widely reported that numerous factors could alter characteristics of gut microbiota at a subject-wise level, e.g., feeding practices, mode of deliver, taking antibiotics and so on, which in turn lead to alterations of subject specific microbiota trajectories. We thus extracted features of longitudinal trajectories of gut microbes at the subject level to fill this knowledge gap. Through our analysis, a comprehensive understanding of the potential associations between characteristics of gut microbiota maturation and development of cognition and temperament during the first three years of life can be gained.

Results

All study activities were approved by the University of North Carolina at Chapel Hill (UNC) and University of Minnesota (UMN) Institutional Review Boards, respectively. Typically developing children meeting the inclusion/exclusion criteria were enrolled (Supporting materials). An accelerated longitudinal design was employed for which the scheduled visits varied among subjects spanning over the first three years of life (Fig. S1). To uncover the potential associations of the developing gut microbiota in the development of cognition/temperament during the first three years of life, the Mullen Scales of Early Learning (MSEL)^32,33 was used to assess cognition for all children and the Infant Behavior Questionnaire Revised (IBQ-R)³³ and Early Childhood Behavior Questionnaire (ECBQ)¹⁹ were used for temperament in children younger and older than 15 months of age, respectively. The MSEL includes five subdomains: gross motor (GM), fine motor (FM), expressive language (EL), receptive language (RL), and visual reception (VR). In addition, an early learning composite (ELC) score was derived for each subject as the summation of FM, EL, RL, and VR. The three composite scales for IBQ-R³³ (Surgency (SUR), Negative Affectivity (NEG), and Orienting/Regulation (REG)) and ECBQ¹⁹ (Negative Affectivity (NEG), Surgency/Extraversion (SUR), and Effortful Control (EFF)) were also employed in our data analyses (Table S1). The demographic information including maternal education, delivery mode, and feeding practices were collected for the parents of the participating infants (Table 1).

Table 1. Demographic information.

Demographic Information	Number	Mean (SD)/Count (% of total)
Infants	N = 374
Infant Age in months		13.0 (8.6)
Male		180 (48.1%)
Female		194 (51.9%)
Delivery mode	N = 366
Vaginal		258 (70.5%)
C-Section		108 (29.5%)
4 Month feeding practice1	N = 271
Exclusively breast fed		190 (70.1%)
Others		81 (29.9%)
Parity	N = 267
1		58 (21.7%)
> 1		209 78.3%)
Maternal education	N = 368
Less than college		49 (13.3%)
College		122 (33.2%)
Graduate school		197 (53.5%)

¹Defined as infants fed less than 4 teaspoons or 20 g per day of non-formula and complementary foods/liquids (water, apple juice, etc.)

Gut microbial composition diversity and cognition/temperament

Fecal samples were longitudinally collected (Figs. 1A, S1) using the OMR-200 kit and analyzed using Next-generation sequencing yielding microbial information at the species level. The temporal behaviors of relative abundance of gut microbiota at the phylum level were consistent with the widely reported temporal maturation of the gut microbiota during infancy (Fig. 1B); the relative abundance of Actinomycetota decreases while the Bacteroidota and Bacillota increase with age¹⁵. We first determined if the Shannon index, Chao1 richness index, evenness, and beta diversity (Fig. 1C–F) were associated with cognition/temperament. Specifically, we regressed out the age effects of the four diversity measures and association analyses were conducted between residues of the four diversity measures and MSEL t-scores and age effect regressed IBQ-R/ECBQ composite scores (Fig. S2). All association analyses discerning the relation between gut microbiota characteristics and cognition/temperament controlled for maternal education and mode of delivery (70.5% born vaginally, Table 1). In addition, since a large majority of infants in our cohort were exclusively breast fed at 4 months of age (70.1%), defined as infants fed less than 4 teaspoons or 20 g per day of non-formula and complementary foods/liquids (water, apple juice, etc.), all statistical analyses controlled for four-month feeding practice. Finally, significant sex differences were observed with MSEL, but not IBQ-R/ECBQ (Tables S2a–S2c). Therefore, sex was further controlled for association analyses between MSEL and microbiota features. The number of subjects and observations for each of the outcome measures are shown in Fig. S3.

Fig. 1 [Images not available. See PDF.]

The distribution of numbers of fecal samples collected at different ages (A). The temporal behaviors of relative abundance of gut microbiota at phylum level (B). The temporal characteristics of gut microbiota diversity measures, including Shannon (C), Chao 1 (D), evenness (E) and beta diversity (F), respectively. A two-piece linear model with subject-wise random intercept was fitted to the trajectories of Shannon (C), Chao1 (D), and beta diversity (F). Using the smallest Akaike information criterion (AIC) as the selection criterion for the transition ages from “rapid” to “slow,” we found that the age transitions were 16 months old for Shannon and 15 months old for both Chao1 and beta diversity. Therefore, we defined 0–15 months as the rapidly changed phase of the three measures (blue dashed lines). We extracted subject-wise random intercept and slope of age effects for all subjects who had at least 2 visits during 0–15 months old using a linear mixed effect model (478 observations from 143 unique subjects). The random slopes of beta diversity in relation to ELC showing that a faster decreased beta diversity (negative slopes) at 12–15 months of age is associated with a higher ELC (G and H).

We found that the Chao1 index was positively associated Early Learning Composite (ELC, Fig. S4, p = 0.022, effect = 0.061, 95% CI (0.007, 0.116)), receptive language (RL, adjusted p = 0.01, effect = 0.069, 95% CI (0.058, 0.080)) and gross motor (GM, adjusted p = 0.045, effect = 0.060, 95% CI (0.019, 0.101)), suggesting a higher richness (alpha diversity) of gut microbiota is associated with a better cognitive ability. However, no associations with temperament were observed after correcting for multiple comparisons.

Our findings appear different from those reported in the literature evaluating associations between gut microbiota and temperament/cognition in infancy. Most of the reported results have been associated with temperament, including an increased Shannon diversity at 1 year old and an elevated risk of behavioral problems at 2 years old (p = 0.087)²⁸, a positive relation between phylogenetic diversity and ECBQ SUR at 18–27 months old²⁷, and a negative association between alpha diversity and negative emotion and fear reactivity at 2.5 months old²⁵. Instead, we observed associations with cognition (ELC, RL and GM). These findings contrast a previously reported inverse relation between different alpha diversity indices measured at 1 year-old and cognitive ability at 2 year-old²⁶ and the lack of association between Shannon diversity and GM and FM, communication and problem-solving skills in another study³⁴. Although several possible explanations could account for the observed discrepancies between our findings and those reported in the literature, three plausible factors are the differences of infants’ ages, study designs and instruments used to assess cognition/temperament. The previously reported studies focused on a relatively short age range of few months^25,27 and/or used a relatively sparse schedule for fecal sample collections^26,28, making it difficult to fully capture how the dynamic maturation phase of both alpha and beta diversities during infancy relate to cognition. In addition, instead of using the MSEL, IBR-R and ECBQ, CBCL and ASQ-3 were employed by Loughman et al. and Sordillo et al., respectively^28,34.

To further confirm our observations, we attempted to determine the age range exhibiting rapid changes of diversity measures. As shown in Fig. 1, the shannon index (Fig. 1C), chao1 index (Fig. 1D), and beta diversity (Fig. 1F) changed rapidly initially and stabilized afterwards. Therefore, we determined the age periods exhibiting “rapid” changes of Shannon, Chao1, and beta diversity using a piecewise linear model and using the smallest AIC as the selection criterion for the transition ages. Please see SI: Determining the rapid and slow phases of diversity measures for additional details. We found that 0–15 months was a rapidly changing time for the three measures (Table S3) dashed blue lines, Fig. 1C,D,F) and thus subject-wise random intercept and slope were extracted during this period. The potential associations between the random slopes of Shannon, Chao1 and beta diversity and cognition/temperament of the entire cohort and within age bins including 0–6, 6–12, 12–15, 15–18, 18–24, and older than 24 months old were determined.

When considering all subjects, no significant relationship was observed. In contrast, for subjects between 12–15 months a faster decrease of beta diversity was associated with higher ELC (p = 0.034, e = -2301.35) (Fig. 1G). Consistent findings were observed by comparing the cognitive/temperament scores between subjects in the top and bottom 25% random slopes; ELC is higher in subjects with faster decrease of beta diversity (difference = 8.10, 95% CI (1.58, 14.62), p = 0.015) for subjects between 12 and 15 months old (Fig. 1H). Results on marginal effects for ELC, Mullen subdomains and different age bins are provided in Figure S5.

Features of gut microbial communities and cognition/temperament

Representative temporal patterns of the eight highest relative abundant species are shown in Fig. 2 (top panel), demonstrating the temporal variabilities of the relative abundance of different species and the substantial variations among subjects (bottom panel). The temporal patterns of the remaining species are shown in Figs. S6–S10. While these temporal patterns appear species specific, gut microbes show strong interdependencies and interactions with either competition or synergies may affect their host’s physiology or vice versa^35,36. In particular, the gut microbiota maturation has been characterized by sequential progression of different community structures during infancy^15,31. Such characteristics cannot be captured fully by microbiota diversity measures. We thus employed network-based approaches where each species was considered as a node and the temporal interactions among pairs of species were represented by edges connecting pairs of species, forming a gut microbial network. A total of 116 species (nodes) with mean relative abundance larger than 0.1% were included to mitigate the need of using zero-inflated approaches and the sum of the relative abundance of the remaining 966 species was included as the ‘Others’ species, forming a microbial network comprising a total of 117 nodes.

Fig. 2 [Images not available. See PDF.]

Temporal characteristics of the relative abundance of the eight highest relative abundance species (upper panel) and the corresponding subject-wise trajectories of these eight species (lower panel) where different colors represent different subjects.

We first computed the correlation matrix for microbiome species using the SparCC³⁷ method as suggested by Weiss et al.³⁸ owing to its ability to infer linear relationships with high precision for high diversity compositions. Subsequently, the modular structures of the microbial community were determined using a classic Louvain algorithm, which parcellated species into non-overlapping modules by maximizing modularity³⁹. A total of six modules were detected for 115 nodes (Fig. 3A), which included 57, 20, 2, 4, 30, and 2 species for modules 1–6, respectively (Table S4). Two microbes were completely isolated from all others and thus excluded. Modules 1 and 5 were dominated by species from the family Lachnospiraceae (Module 1: 21.5% and Module 5: 26.7%), whereas module 2 contained primarily Bacteroidaceae (55%) (Fig. S11). Further examining the members of Modules 1 and 5 at the genus level, they shared Blautia, Lachnospiraceae_u_g and Roseburia, while Anaerostipes, Coprococcus, and Lachnoclostridium were unique to Module 1 and Dorea and Tyzzerella for Module 5, respectively (Table S4).

Fig. 3 [Images not available. See PDF.]

The gut microbial network structure where each circle represents a species, colors represent different modules, and the lines indicate the connections among species (A). The relative abundance of species with age for each module are shown in (B) for modules 1–6, respectively. The mean interaction strength with age for each module are shown in (C) for modules 1–6, respectively. Finally, the heatmap showing the associations between the mean relative abundance and interaction strength with cognition and temperament is provided in (D). The symbols in D indicate adjusted p-values: “!”: p < 0.001; “$”: 0.005 <  = p < 0.01; “*”: 0.01 <  = p < 0.05 and the color bar indicates effect sizes.

We then analyzed the temporal behaviors of the mean relative abundance (Fig. 3B) and interaction strength⁴⁰ (Figs. 3C) of these modules and examined how they may be associated with cognition/temperament. Mean relative abundance of module 1 was negatively associated with ELC (effect = − 6.18, adjusted p = 0.0176, 95% CI [− 11.27, − 1.09]) and RL (effect = − 8.18, adjusted p = 0.000486, 95% CI [− 12.27, − 4.08]) (Fig. 3D) whereas module 2 (Bacteroidaceae) showed a positive association with ELC (effect = 7.94, adjusted p = 0.0174, 95% CI [1.40, 14.48]). In contrast, a greater number of associations, encompassing both cognition and temperament, were observed when examining the interaction strength within each module. The species interaction strengths of Modules 3 (ELC: effect = − 1.79, adjusted p = 0.0060, 95% CI [− 3.06, − 0.52]) and FM: effect = − 1.61, adjusted p = 0.0103, 95% CI [− 2.63, − 0.59]) and 5 (GM: effect = − 12.16, adjusted p = 0.0443, 95% CI [− 21.24, − 3.08] were associated with cognition while Modules 1 (ECBQ_EFF: effect = − 5.69, adjusted p-value = 0.0257, 95% CI [− 9.84, − 1.54]) and 4 (ECBQ_NEG: effect = 0.56, adjusted p = 0.0205, 95% CI [0.16, 0.95]) with temperament (Fig. 3D). Note that the reported effect sizes above as well as hereafter are raw regression coefficients (β) using the linear mixed-effects models. These represent the change in cognitive/temperament scores per unit change in the microbial feature. We have included effect size (beta/standard deviation) in the Supplementary Table S7, allowing direct comparison of the magnitude of associations across variables (Table 2).

Table 2. The 10 pairs of species exhibit the highest positive (co-occurrence) and negative (co-exclusion) connections.

Number	Co-occurrence		Co-exclusion
	Specie 1	Specie 2	Specie 1	Specie 2
1	Collinsella sp. 4_8_47FAA	Collinsella aerofaciens	Ruminococcus sp. 5_1_39BFAA	Ruminococcus sp. SR1/5
2	Anaerostipes hadrus	Lachnospiraceae bacterium 5_1_63FAA	Blautia sp. KLE 1732	Others
3	Clostridiales bacterium VE202-03	Lachnospiraceae bacterium 7_1_58FAA	Roseburia inulinivorans	Others
4	Anaerostipes hadrus	Clostridium sp. SS2/1	Lachnospiraceae bacterium 5_1_63FAA	Others
5	Veillonella parvula	Veillonella dispar	butyrate producing bacterium SS3/4	Others
6	Clostridiales bacterium VE202-18	Erysipelatoclostridium ramosum	Anaerostipes hadrus	Others
7	Lachnospiraceae bacterium 5_1_63FAA	Clostridium sp. SS2/1	Bifidobacterium longum	Roseburia inulinivorans
8	Flavonifractor plautii	Lachnospiraceae bacterium 7_1_58FAA	Eubacterium hallii	Others
9	Erysipelotrichaceae bacterium 6_1_45	Erysipelotrichaceae bacterium 2_2_44A	Blautia obeum	Others
10	Flavonifractor plautii	Clostridiales bacterium VE202-03	Ruminococcus torques	Others

Furthermore, we identified key species within the microbial network by calculating the degree of each node (species), where the degree of a node denotes the total number of connections it maintains with other species (Fig. 4A). That is, a species possesses a high degree would be indicative of a highly interconnected species in the network. To evaluate the impact of removing highly connected species on the network’s interaction efficiency, we performed simulations of targeted attacks—sequentially removing the highest-degree species, and random attacks—removing species indiscriminately. After each removal, we recalculated the network’s global efficiency, which provided insights into the overall interaction and functional integration among the microbes in the human gut. Our results revealed a significant reduction in global efficiency following targeted attacks (removing nodes sequentially based on their degrees in a descending order) compared to random attacks (Fig. 4B). This finding indicates that the removal of highly connected species led to a rapid reduction in global efficiency, suggesting that the microbes with a high degree play key roles for facilitating functional interaction among gut microbes.

Fig. 4 [Images not available. See PDF.]

The degrees for all species (A) where the red lines indicate the 25 highest degrees. The effects of simulated targeted and random attacked on the global efficiency of the microbial network (B). The associations between MSEL, IBQ-R, and ECBQ with the Microbiome network features are shown in (C). The symbols in C indicate adjusted p-values: “!”: p < 0.001; “#”: 0.001 <  = p < 0.005; “$”: 0.005 <  = p < 0.01; “*”: 0.01 <  = p < 0.05 and the color bar indicates effective sizes). The x axis in C indicate the 25 degree hubs, 10 strongest co-occurrence pairs, and 10 strongest co-exclusion pairs.

To determine how high degree species are associated with cognition/temperament, we focused our analyses on the 25 highest degree species, defined as the degree hubs (degree cutoff as 53, red lines, Fig. 4A). Of these 25 hubs (Table 3), most of the hubs were members of the Module 1 (15) and Module 5 (8), and 2 hubs in Module 4 (Fig. 3A). Since the “Others” node containing 966 species inevitably exhibited the highest degree, was excluded from the subsequent association analysis. The remaining hubs are in the order Clostridiales (14), Bifidobacteriales (4), Lactobacillales (2), Veillonellales (2), Enterobacterales (1), and Erysipelotrichales (1) (Table 3). Finally, the B longum exhibited the highest degree (69), suggesting that it is highly connected with other microbes (Table 3).

Table 3. The 25 degree hubs.

Kindom	Phylum	Class	Order	Family	Genus	Species	Degree	Module
Others							76	1
Bacteria	Actinobacteria	Actinobacteria	Bifidobacteriales	Bifidobacteriaceae	Bifidobacterium	Bifidobacterium longum	69	5
Bacteria	Proteobacteria	Gammaproteobacteria	Enterobacterales	Enterobacteriaceae	Escherichia	Escherichia coli	65	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Blautia	Ruminococcus torques	65	5
Bacteria	Firmicutes	Bacilli	Lactobacillales	Enterococcaceae	Enterococcus	Enterococcus faecalis	64	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Blautia	Blautia sp. KLE 1732	63	5
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Blautia	Blautia obeum	63	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Roseburia	Roseburia inulinivorans	61	5
Bacteria	Firmicutes	Clostridia	Clostridiales	Eubacteriaceae	Eubacterium	Eubacterium hallii	60	1
Bacteria	Firmicutes	Negativicutes	Veillonellales	Veillonellaceae	Veillonella	Veillonella dispar	60	1
Bacteria	Actinobacteria	Actinobacteria	Bifidobacteriales	Bifidobacteriaceae	Bifidobacterium	Bifidobacterium bifidum	59	1
Bacteria	Firmicutes	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus	Streptococcus sp. C150	59	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Anaerostipes	Anaerostipes hadrus	58	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Coprococcus	Coprococcus comes	58	1
Bacteria	Actinobacteria	Actinobacteria	Bifidobacteriales	Bifidobacteriaceae	Bifidobacterium	Bifidobacterium breve	57	4
Bacteria	Actinobacteria	Actinobacteria	Bifidobacteriales	Bifidobacteriaceae	Bifidobacterium	Bifidobacterium dentium	57	4
Bacteria	Firmicutes	Clostridia	Clostridiales	Ruminococcaceae	Ruminococcus	Ruminococcus lactaris	57	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Lachnospiraceae_u_g	Lachnospiraceae bacterium 5_1_63FAA	55	5
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Roseburia	Roseburia hominis	55	1
Bacteria	Firmicutes	Negativicutes	Veillonellales	Veillonellaceae	Veillonella	Veillonella parvula	54	5
Bacteria	Firmicutes	Clostridia	Clostridiales	Clostridiales_u_f	Clostridiales_u_g	butyrate producing bacterium SS3/4	54	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Lachnospiraceae	Dorea	Dorea longicatena	54	5
Bacteria	Firmicutes	Clostridia	Clostridiales	Oscillospiraceae	Oscillospiraceae_u_g	Oscillospiraceae bacterium VE202-24	54	5
Bacteria	Firmicutes	Erysipelotrichia	Erysipelotrichales	Erysipelotrichaceae	Erysipelatoclostridium	Erysipelatoclostridium ramosum	54	1
Bacteria	Firmicutes	Clostridia	Clostridiales	Ruminococcaceae	Ruminococcus	Ruminococcus bromii	53	1

We observed more associations between the degree hubs and cognition than that with the temperament. In particular, RL (Fig. 4C) was associated with multiple degree hubs, including a negative association with B. longum (effect = − 7.41, adjusted p = 0.000195, 95% CI [− 10.92, − 3.90]), V. dispar (effect = − 190.44, adjusted p = 0.0117, 95% CI [− 312.27, -68.62]) and V. parvula (effect = − 139.41, adjusted p = 0.0384, 95% CI [− 241.57, − 37.25]), but was positively associated with Butyrate-producing bacterium SS3/4 (effect = 393.74, adjusted p = 0.00109, 95% CI [186.27, 601.20]), E. hallii (effect = 162.96, adjusted p = 0.0290, 95% CI [47.17, 278.76]), R. bromii (effect = 91.43, adjusted p = 0.00110, 95% CI [43.11, 139.75]), Roseburia hominis (effect = 206.40, adjusted p = 0.0433, 95% CI [52.66, 360.15]), Roseburia inulinivorans (effect = 149.26, adjusted p = 0.0300, 95% CI [43.92, 254.59]), and R. lactaris (effect = 291.08, adjusted p = 0.0241, 95% CI [88.35, 493.82]). In addition, ELC was positively associated with R. bromii and R. lactaris (effect = 63.30, adjusted p = 0.0387, 95% CI [3.01, 123.58]; effect = 335.56, adjusted p = 0.00805, 95% CI [88.56, 582.57]) whereas V. dispar (effect = − 202.54, adjusted p = 0.00833, 95% CI [− 352.36, − 52.73]) and V. parvula (effect = − 156.86, adjusted p = 0.0144, 95% CI [− 283.98, − 29.73]) were negative. Finally, it is worth noting that with the exception of EL, R. bromii also exhibited significant associations with all remaining subdomain scores, positive with RL and VR (effect = 65.93, adjusted p = 0.0350, 95% CI [18.11, 113.76]) and negative with FM (effect = − 65.24, adjusted p = 0.0343, 95% CI [− 112.19, − 18.29]) and GM (effect = − 69.39, adjusted p = 0.0322, 95% CI [− 119.14, − 19.64]).

In contrast to the numerous associations between degree hubs and cognitions, only two associations with temperament were observed; B. KLE 1732 was positively associated with IBQR_REG (effect = 11.81, adjusted p = 0.0364, 95% CI [2.69, 20.94]) and the meanhub with the ECBQ_NEG (effect = 43.62, adjusted p = 0.0187, 95% CI [13.09, 74.15]). It should be noted that the aforementioned associations between microbes and cognitions are largely for species in the order Clostridiales and Veillonellales, respectively.

Most species within module 2 and module 5 are positively correlated (59/67 and 163/195 respectively), suggesting that the relative abundance of these species increase or decrease together (co-occurrence)³⁶. Conversely, module 1 exhibits a notable proportion of negative correlations, with about 48% (544/1139) of interactions showing negative associations, suggesting competing or inhibiting relations among these species (co-exclusion)³⁶. Modules 4 and 6 show predominantly positive correlations, suggesting a cooperative or mutualistic relationship. Finally, Module 3 contains two microbes (Enterococcus faecium and Anaerostipes sp. 3_2_56FAA, Table S4), which were negatively interacted, indicating an inhibitory interaction between the two species. To determine if the interaction patterns among pairs of species may be associated with cognition/temperament, we identified the 10 most positively and negatively correlated species-pairs (Table 2), respectively. The sums/differences of the relative abundance of each of the 10 co-occurrence/co-exclusion³⁶ pairs were calculated. Specifically, for positively correlated pairs (co-occurrence), the sums of relative abundances were computed to reflect their combined abundance. For negatively correlated pairs (co-exclusion), the absolute differences in relative abundances were calculated to represent the extent of their exclusivity. Interestingly, unlike the relative abundance of the identified degree hubs which exhibited more associations with cognition than that of temperament, the interaction strengths were more associated with temperament. We found (Fig. 4C) that the co-occurrence of Flavonifractor plautii—C. bacterium VE202-03 (pos_10) and Flavonifractor plautii—Lachnospiraceae bacterium 7_1_58FAA (pos_8) were positively associated with ECBQ NEG (effect = 52.94, adjusted p = 0.0325, 95% CI [12.05, 93.83]; effect = 64.75, adjusted p = 0.00938, 95% CI [22.74, 106.76]) and negatively associated with IBQR NEG (effect = − 16.88, adjusted p = 0.0323, 95% CI [− 29.63, − 4.14]; effect = − 19.07, adjusted p = 0.0171, 95% CI [− 32.31, − 5.84]). Since IBQ-R and ECBQ were collected for children younger and older than 15 months of age, respectively, these findings suggest that their associations with the NEG trait are consistent throughout the first three years of life. In contrast, co-exclusion of R. sp. 5_1_39BFAA—R. sp. SR1/5 (neg_1) was negatively associated with IBQR-REG (effect = − 4.22, adjusted p = 0.0259, 95% CI [− 7.33, − 1.10]) and B. sp. KLE 1732—“Others” (neg_2) was positively associated with IBQR SUR (effect = 2.02, adjusted p = 0.0307, 95% CI [0.50, 3.54]).

Regarding cognition, co-occurrence of V. parvula—V. dispar (pos_5) was negatively associated with ELC (effect = − 93.49, adjusted p = 0.00918, 95% CI [− 163.40, − 23.58]) and RL (effect = − 85.24, adjusted p = 0.0178, 95% CI [− 142.82, − 27.65]). Finally, co-exclusion of B. longum and R. inulinivorans (neg_7) was negatively associated with RL (effect = − 7.45, adjusted p = 0.000150, 95% CI [− 10.93, − 3.97]).

Three main findings using network-based approaches warrant additional discussion (Table 4). First, similar to diversity measures, the network-derived features of the gut microbiome showed more frequent and robust associations with cognitive outcomes than with temperament measures, although some associations with temperament were also detected. For instance, at the module level, the mean relative abundance of Module 1 showed negative associations with RL and ELC, while Module 2 exhibited a positive association with ELC, collectively highlighting that different modules and species may contribute distinctly to cognitive domains. These findings underscore that beyond simple diversity metrics, the network structure and species interactions of the gut microbiome may be influential in cognitive development.

Table 4. Summary of the main findings of network-based analysis.

			MSEL						IBQ-R			ECBQ
			ELC	GM	FM	EL	RL	VR	SUR	NEG	REG	SUR	NEG	EFF
Module	Relative Abundance	Module 1	−				−
		Module 2	+
	Interaction Strength	Module 1												−
		Module 3	−		−
		Module 4											+
		Module 5		−
Mean Hub							−						+
Degree Hub	BPB SS3/4						+
	B. KLE 1732										+
	E. hallii						+
	OB VE202−24				−
	R.bromii		+	−	−		+	+
	R.hominis						+
	R.inulinivorans						+
	R.lactaris		+				+
	V.dispar		−				−
	V.parvula		−				−
Co-occurrence	pos_5 (V.parvula and V.dispar)		−				−
	pos_8 (F.plautii and L.bacterium 7_1_58FAA)									−			+
	pos_10 (F.plautii and C.bacterium VE202−03)									−			+
Co-exclusion	neg_1 (R.sp. 5_1_39BFAA and R.sp. SR1/5)										−
	neg_2 (B.sp. KLE 1732 and Others)								+
	neg_7 (R.torques and Others)						−

Second, our results revealed complex patterns of species co-occurrence and co-exclusion that were linked to cognitive and temperament metrics. While Modules 2 and 5 featured predominantly positive species-to-species correlations—indicative of cooperative or synergistic interactions—Module 1 contained a notable proportion of negative correlations, suggesting competition or inhibitory interactions among its members. Strong co-occurrence relationships, such as the V. parvula–V. dispar pair (pos_5, Fig. 4C), were negatively associated with MSEL ELC and RL. Furthermore, both V.dispar and V. parvula were among the top 25 degree species (degree hubs), which were also negatively associated with ELC and RL (Fig. 4C). These findings suggest that V.dispar and V. parvula not only individually but their interaction were negatively associated with ELC and RL. Nevertheless, these associations appear differ from results reported in the literature. Decreased relative abundance of the genus Veillonella was previously observed in autistic children⁴¹ and in PD with mild cognitive impairment⁴². These observed discrepancies are most likely attributed by the lack of results focusing on the species level. Our findings underscoring the importance to evaluate these potential associations with cognition and temperament at the species level. Finally, our findings lend support to the idea that considering the microbial community as an interconnected network can unveil associations that remain hidden when focusing solely on single species or standard diversity indices.

Third, examining the network’s topology and identifying key species hubs shed light on the functional importance of highly connected taxa. We found that removal of high-degree species led to a substantial decrease in global efficiency, emphasizing their central roles in sustaining network integrity. The top 25 hubs—which included B. longum (the highest degree hub) and several Ruminococcus, Veillonella, and Enterobacterales species—showed multiple associations with cognition and temperament. For example, B. longum, a commonly recognized beneficial gut bacterium, along with other hub species in Clostridiales and Veillonellales, was linked to RL as well as other cognitive domains. The negative associations of pos_5 (V. parvula–V. dispar co-occurrence) and neg_1 (Ruminococcus sp. pairs) with specific cognitive measures highlight the dynamic and nuanced influence of microbial relationships during infancy and early childhood. Shifts in microbial composition and interaction patterns over time may contribute differently to cognitive or temperamental outcomes depending on a child’s developmental stage.

Fourth, one of our major findings is the associations between R. bromii and cognition, including positive associations with comp, RL and VR and negative with FM and GM, respectively. R. bromii is known to promote the activity of butyrate-producing bacteria ⁴³ and produce mostly acetate as SCFAs, which in turn could enhance intestinal barrier function and reducing inflammatory conditions. In addition, Ze et al⁴⁴ reported that R. bromii plays a pivotal role for the degradation of resistant starch in the human large intestine, which in turn could potentially provide health benefits on reducing insulin resistance and preventing colorectal cancer. Nevertheless, the negative associations of R. bromii with motor functions (FM and GM) remain unclear and warrant additional studies.

Finally, the relative abundance of B. longum as well as its interaction with R. inulinivoranis (co-exclusion) were negatively associated with RL. These findings are puzzling and contradicts to the consensus that B. longum is a bacterium associated with healthy conditions. Although many factors could potentially explain the negative findings in our study, it should be noted that most of the interventional studies reporting cognitive benefits of B. longum largely focused on adults with relatively few studies on early infancy as reviewed by Eastwood et al⁴⁵. Furthermore, a combined supplement of probiotics including B. longum from 12 to 24 months old showed no cognitive benefits assessed at 24 months⁴⁶. Finally, to the best of our knowledge, our findings are the first reported results on the associations between B. longum and language development during early infancy.

Features reflecting subject-wise gut microbial trajectories and cognition/temperament

One of the unique features of our study is the longitudinally collected fecal samples 0–3 years old (Figure S1, Table S5), enabling evaluation of how subject-wise microbiota trajectories could be associated with cognition/temperament. In particular, the developmental trajectories not only vary among microbial species but also between children (lower panel Fig. 2, Figures S7–S10). Nevertheless, general temporal trajectory patterns emerge; some species exhibit a low or high relative abundance while some exhibit a low relative abundance during infancy and transition to a high relative abundance in toddlerhood and vice versa (Fig. 2). A Gaussian mixture model with three subject-wise random effects, high- (HA) and low-abundance (LA) groups and one logistic model for the probability of a sample belonging to high abundance (PHA) group was used to fit the trajectory of each species. To ensure more accurate estimates of the random effects, only subjects with at least three visits (Figure S12) and species with a mean relative abundance > 0.5% across all subjects (31 species) were included for subsequent analyses. A representative example of the subject-wise relative abundance trajectories of Bacteroides_u_s and the corresponding three random effects are shown in Figs. 5A and 5B and the remaining species are shown in Figures S13–S17. A hierarchical clustering approach was used based on the above three subject-wise random effects and the average silhouette distance was used to select the optimal number of clusters for each species. Four clusters were identified for Bacteroides_u_s (Figs. 5C–5F) and significant differences of ECBQ motor activation (p = 0.009, Fig. 5G) and sadness (p = 0.002, Fig. 5H) were observed among subjects in different clusters. Figure 5I summarizes the clustering results for all species. Among the species that exhibited significant differences between clusters, species in the Ruminococcaceae family exhibit significant differences in cognition and Bacteroidaceae and Lachnospiraceae show differences in ECBQ subdomain scales. In contrast, species show cluster differences for the IBQ-R span in many of the families. The most significant differences (raw p < 0.005) include four ECBQ and one MSEL scales. In addition, Bifidobacterium adolescentis and unclassified Ruminococcus sp. 5_1_39BFAA were significantly different for ECBQ Activity (p = 0.0006) and Soothability (p = 0.008), respectively. Furthermore, the clusters of Bacteroides vulgatus exhibit significant differences of MSEL visual reception (p = 0.009). Note that most of the associations were with temperament traits instead of cognition, contrary to the results using diversity and network-based microbiota features as described above.

Fig. 5 [Images not available. See PDF.]

The subject-wise relative abundance trajectories of Bacteroides_u_s (A) and the three random effects (B) where brown, blue and green lines represent high abundance (HA), low abundance (LA) and probability of transition to high abundance (P_HA), respectively. The clustering results of the Bacteroides_u_s where four clusters (C-F) were obtained representing subjects with different patterns of trajectories. Significant differences of ECBQ motor activation (G) and sadness (H) were observed among subjects in the four clusters. The associations between MSEL, IBQ-R, and ECBQ and the cluster results are shown in (I). Different shapes represent the statistical results for clustering analysis: squares: p < 0.001; circles: 0.001 < p < 0.005; triangles: 0.005 < p < 0.01; diamonds: 0.01 < p < 0.05.

We further evaluated how each of the three subject-wise random effects was associated with cognition/temperament (Figs. 6A-C). Evidently, the associations were feature dependent; more associations were observed with the PHA and HA than that of LA features (Table 5). In addition, consistent with the results using the clustering approach, there were more and stronger associations with temperament (IBQ-R and ECBQ) than that of cognition (Table 5).

Fig. 6 [Images not available. See PDF.]

The associations between MSEL, IBQ-R, and ECBQ and with three subject-wise random effects, including P_HA (A), HA (B) and LA (C), respectively. Different symbols indicate the statistical results for the association analyses with the three subject-wise random effects: “!”: p < 0.001; “#”: 0.001 < p < 0.005; “$”: 0.005 < p < 0.01; “*”: 0.01 < p < 0.05. The color bar indicates log(p-values).

Table 5. Number of associations with cognition/temperament at different statistical significance levels of three trajectory features.

		MSEL	IBQ-R	ECBQ	Total
High abundance (HA)	P < 0.001
	0.001 < P < 0.005		9	1	10
	0.005 < p < 0.01	1	6	1	8
	0.01 < p < 0.05	4	23	33	60
Low abundance (LA)	P < 0.001		1	1	2
	0.001 < P < 0.005		2	4	6
	0.005 < p < 0.01		5	5	10
	0.01 < p < 0.05	5	20	26	51
Probability of Transition (PHA)	P < 0.001			1	1
	0.001 < P < 0.005	1	1	3	5
	0.005 < p < 0.01		3	4	7
	0.01 < p < 0.05	4	18	29	51

Regarding how the three random effects were associated with cognition/temperament, the strongest associations among the three features were in genus Bacteroides and with ECBQ, including, HA: Bacteroides_u_s and ECBQ Motor activation (p = 0.001, +); LA: Bacteroides ovatus and ECBQ Soothability (p = 0.0002, −); and PHA: Bacteroides dorei and ECBQ Soothability (p = 0.0008, +) where the sign in the parentheses indicate the direction of association. The finding of Bacteroides_u_s is consistent with that using the cluster approach (Fig. 5G). Furthermore, both Bacteroides ovatus and Bacteroides dorei are associated with the ECBQ Soothability, suggesting that subjects > 15 months old whose gut microbiota containing Bacteroides dorei with PHA of transitioning to HA group as well as subjects with a lower Bacteroides ovatus abundance in the LA group were associated with a higher ECBQ Soothability. Finally, although the associations were not as strong as those with temperament, it is worth highlighting the two strong associations with the RL, including unclassified Ruminococcus sp. 5_1_39BFAA (HA: p = 0.0056, +) and Blautia sp. KLE 1732 (PHA, p = 0.002, −). These findings suggest that the subjects who exhibited a high abundance of Ruminococcus sp. 5_1_39BFAA in the HA group were associated with higher RL scores. Furthermore, subjects who exhibited a higher probability of belonging to high abundance group for Blautia sp. KLE 1732 were associated a lower RL scores.

Finally, several species in genus Bifidobacterium also exhibited strong associations with temperament traits. Bifidobacterium is one of the most well studied bacteria and widely implicated to possess health benefits. Our findings are largely consistent with the reported positive health benefits of Bifidobacterium but provide new insights into the differential associations with temperament traits among the species in genus Bifidobacterium. Major associations (raw p < 0.005) between temperament traits with B. bifidum, B. longum, B. pseudocatenulatum, and B. adolescentis were observed (Table S6). The PHA of B. bifidum exhibited an inverse association with ECBQ Sociability (p = 0.0022) and Attentional shifting (p = 0.0048), implying that children > 15 months old and exhibiting a high probability of transition to high relative abundance were associated with lower ECBQ Sociability (enjoying interaction with others) and Attentional shift (the ability to shift attention from one task to another) scores. In addition, children < 15 months old, in the HA group and exhibiting a higher relative abundance of B. bifidum were inversely associated with IBQ-R Sadness whereas B. longum were positively with Vocal reactivity: vocalization exhibited by infants. As discussed above, the relative abundance of both B. bifidum and B. longum is higher than other species until about 15 months and markedly reduced > 15 months (Figure S6). Similar findings are also observed for B. pseudocatenulatum, and B. adolescentis. Our findings further confirm the importance of age-appropriate gut microbiota for temperament traits.

Discussion and conclusion

Here we provided evidence on the possible links between the development of gut microbiota and cognition and temperament traits in typically developing children 0—3 years old, an age period of rapid maturation of both biological systems. Our findings not only illustrate that associations between gut microbiota and temperament/cognition varied with the analytical approaches but also highlighted differential gut microbial features in association with cognition and temperament traits. Furthermore, our findings signify the importance of age-appropriate gut microbiota in associations with temperament traits.

The causal link of the gut microbiota with brain related outcomes is best illustrated in mice models. Comparisons between germ-free (GF) and specific pathogen-free (SPF) mice show a higher plasma corticosterone and adrenocorticotropic hormone ⁴⁷ in GF mice than that of SPF mice under stress. In addition, increased motor activity, reduced anxiety ⁴⁸, and significant social impairment ⁴⁹ were observed in GF over that of SPF mice. Unlike a strong body of preclinical evidence elucidating the links between gut microbiota and social behavior, human evidence is relatively limited. Alterations of gut microbiota diversity has been widely observed in autism spectrum disorder (ASD) patients and have been implicated to potentially serve ²² as a biomarker of ASD ^{8, 9, 10–11}. Consistent with these reported results, we observed strong associations in several temperament traits, particularly Soothability and Sociability, further suggesting potential links between gut microbiota and social behaviors. Although we did not directly study social behavior in our cohort, our findings were largely in the language subdomains like RL, which is the pivotal cognitive abilities for communication that could benefit social interaction. Su et al. ⁵⁰ studied the potential links between language functions and social motivation in ASD children 14–31 months old and reported that children with a better language use 2 years after the initial assessments were associated with stronger social motivation. Mulvey and Jenkins ⁵¹ reported that although language alone does not predict social behavior, a strong positive association between the two has been observed in preschoolers. Collectively, our findings reveal novel insights into the potential roles of the gut microbiota in relation to temperament and cognition and converge on the potential interplay of language ability and social behavior during a period when both brain functions and microbiota undergo rapid maturation.

Our study has several limitations. First, although the use of an accelerated longitudinal design offered several advantages, it also hampered our ability to conduct prediction analyses discerning the effects of gut microbiota on cognition and temperament. In addition, since the ages at which subjects were enrolled and the age range of each subject within the study varied across our cohort, associations results for specific age ranges can not be accomplished. Future studies using a classic longitudinal design is warranted to provide insights into age specific associations between gut microbiota and cognitions/temperament. Second, more than 50% of the mothers in our cohort had graduate degrees. Although the potential effects of maternal education were controlled, it is important to interpret our findings by considering the relatively high maternal education levels in our cohort. Third, while the use of whole genome sequencing offers a deeper read of microbiota into the species level, many species have an extremely low abundance level. Therefore, the inference based on these species might be uninformative. To address this issue, we chose a threshold of mean abundance of 0.5% to maintain sufficient information for statistical inference. Fourth, uncorrected associations were reported for subject-wise trajectory analyses owing to the number of microbiome species together with different domains of cognitive and temperament outcomes. Finally, it should be noted that other covariates including diet, both pre- and post-natal environment exposures, antibiotic use, illness, and attendance to day care are potential confounders, and were not controlled in our study. Future work using a classic longitudinal design, a cohort representing the general population and controlling critical covariates will be needed, which calls for a multicenter study capable of mitigating the above limitations. Nevertheless, it should be noted that our study has major strengths, including a relatively large sample size when compared to results reported in the literature, covering a critical age period when our brain and gut microbiota undergo rapid and dynamic maturation, a longitudinal design, and the use of novel statistical approaches to discern associations between gut microbiota and temperament/cognition. Future studies focusing on elucidating the potential causal roles of identified microbes are warranted.

Methods

Subjects were enrolled at both University of North Carolina at Chapel Hill (UNC) and University of Minnesota (UMN) using site-based research registries for identifying research participants between birth and 3 years of age. All research activities were approved by the Institutional Review Boards of UNC and UMN. Written informed consent was obtained from the parents for the participation of both themselves and their infants prior to any study activities. All research was performed in accordance with relevant guidelines/regulations. Details of the inclusion and exclusion criteria and the visiting schedule are provided in the supporting materials.

Stool sample collection and processing

Parents were instructed to collect fecal samples from their child’s diaper at home using the Omnigene Gut sample collection kit (DNA GenoTek, Ontario, Canada) 24 h prior to the in-person visit and bring the collected samples to the in-person visits. If parents forgot to collect fecal samples prior to the scheduled visit, fecal samples were collected on site if a child had a bowel movement during the study appointment. If a child did not have a bowel movement during the scheduled visit, a pre-addressed return envelope was provided and parents were asked to collect the sample from their child’s diaper at home.

Processing the fecal samples

While the collected fecal samples should remain stable in the collection tube for up to 60 days, all collected samples were processed within a week. Specifically, fecal samples were first placed in a dry bead bath to loosen the fecal sample. A sterile transfer pipette was used to transfer the fecal sample into the Eppendorf tubes. The fecal sample was split evenly between 2 × 1.5 ml Eppendorf tubes and was frozen and place in the -80 freezer immediately. Finally, all collected fecal samples were shipped to CosmosID Inc. (Germantown, MD, USA) for further analyses detailed below.

DNA extraction

According to the manufacturer’s protocol, DNA from samples was isolated using the QIAGEN DNeasy PowerSoil Pro Kit (Qiagen, Germantown, MD, USA). Extracted DNA samples were quantified using Qubit 4 fluorometer and Qubit™ dsDNA HS Assay Kit (Thermofisher Scientific, MA, USA).

Library preparation and sequencing

DNA libraries were prepared using the Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA, USA) and IDT Unique Dual Indexes with total DNA input of 1 ng. Genomic DNA was fragmented using a proportional amount of Illumina Nextera XT fragmentation enzyme. Unique dual indexes were added to each sample followed by 12 cycles of PCR to construct libraries. DNA libraries were purified using AMpure magnetic Beads (Beckman Coulter, Brea, CA, USA) and eluted in QIAGEN EB buffer. DNA libraries were quantified using Qubit 4 fluorometer and Qubit™ dsDNA HS Assay Kit. Libraries were then sequenced on an Illumina NovaSeq 6000 System with S4 Flow Cell.

Bioinformatics analysis

Unassembled sequencing reads were directly analyzed by CosmosID-HUB Microbiome Platform (CosmosID Inc., Germantown, MD) described previously ⁵². Briefly, the platform employed curated genome databases together with a high-performance data-mining algorithm that allowed to rapidly disambiguate hundreds of millions of metagenomic sequence reads into the discrete microorganisms engendering the particular sequences.

In terms of the quality checks, upon data generation, raw data were backed up to Amazon AWS and run through fastqc. A multiqc report was generated. The multiqc report was checked to ensure read depth thresholds were met, and that there were no abnormalities with read quality, duplication rates, or adapter content. Taxonomic results were checked on the http://app.cosmosid.com platform to ensure there were no contamination or barcoding issues. The filtering threshold which determines if results are considered significant was based on statistical scores determined by analyzing a large number of diverse metagenomes.

Assessments of cognition and temperament

The Mullen Scales of Early Learning (MSEL)³² were employed to assess cognition for all subjects. In contrast, infant temperament was assessed using the Infant Behavior Questionnaire Revised (IBQ-R)³³ for infants younger than 15 months old while the Early Childhood Behavior Questionnaire (ECBQ)¹⁹ was used for all remaining subjects. More detailed information of the MSEL, IBQ-R, and ECBQ is provided in the Supporting materials.

Microbial composition diversity measures

The calculation of diversity was conducted using the “microbiome” in R package. Specifically, Shannon diversity was calculated as where P_i is the relative abundance of each species. The Chao1 was calculated as where is the observed number of species, and are the numbers of OTUs with only 1 and 2 sequences, respectively. The Chao1 index estimates species richness, accounting for undetected taxa by incorporating singleton and doubleton counts. Finally, beta diversity was calculated as the dissimilarity of each sample against the group mean where the values were then compared between groups to compare the differences in group homogeneity. The correlations among the four parameters was provided in Fig. S18.

Network analysis

Microbiota community can be represented by a network which compromises of nodes and edges, which connect between pairs of nodes. For a weighted network, the edges are weighted according to how strong the connection/interaction between a given pair of nodes. In the context of a microbial network, each node represents an individual microbe (at the species level for our study) and each edge represents the extent of interaction between a given pair of microbes. In this study, network analyses were applied to identify the characteristics of the microbial community during the first 3 years of life. A total of 116 species with mean relative abundance larger than 0.1% and the sum of the relative abundance of the remaining 966 species was included as the ‘Others’ node, leading to a total of 117 nodes for the gut microbial network.

Following the suggestions of Weiss et al.³⁸, the SparCC³⁷ method was used to compute the correlation matrix for microbiome species owing to its ability to infer linear relationships with high precision for high diversity compositions. We utilized the software FastSpar⁵³, which is a C +  + implementation of the SparCC algorithm to efficiently computation of the correlation matrix. We bootstrapped the species abundance dataset by 1000 times and used them as well as the original species to compute the correlation matrix with a correlation strength exclusion threshold of 0.1, a maximum exclusion iterations of 10, and 100 iterations. We then used the 1000 permuted correlation matrixes generated from the bootstrapped dataset to calculate the permuted p-values for the original correlation matrix.

The Louvain algorithm³⁹ (igraph R package function ‘cluster_louvain’) was employed to detect the modular structure of the microbial network with positive and negative weights, using correlations with an absolute value larger than 0.2 and permuted p-value smaller than 0.05. By optimizing the modularity metric, the Louvain algorithm identifies groups of nodes that exhibit denser connections among themselves compared to the rest of the network. Specifically, nodes were assigned to modules by a greedy algorithm that maximizes the modularity index Q, defined as:

1

Here, represents the weight between nodes and , and represent the positive and negative elements of . In addition, , , , , and is the community to which node is assigned, the function is 1 if and 0 otherwise. With the classic Louvain approach, each species was assigned to a module.

Identifying degree hubs

One of the network-based approaches is to identify nodes playing pivotal roles on communicating information throughout the networks. The critical nodes are also known as “hubs” of a given network. Here, we calculated the “degree” of each node based on the binarized adjacent matrix, which represents the sum of the number of edges connected to a given node. Subsequently, all of the nodes in the network were ranked based on the degree of each node. We considered the species within the top 20th percentile of the degree distribution as the degree hubs in our study, leading to a total of 25 degree hubs. This ensures that the identified hubs are those with relatively higher connectivity compared to the overall network.

Determining subject-wise interaction strengths among pairs of nodes within each module

As outlined above, the microbial network was constructed by calculating the correlation between pairwise species from all subjects. Therefore, the interaction strengths of edges represented the group level interaction patterns but not at a subject-wise level. To this end, instantaneous co-fluctuation magnitude between pairs of species were calculated to represent the interactions between two species for each subject⁴⁰. Specifically, the relative abundance level of each species across subjects was z-scored first. Next, for all pairs of species, we calculated the dot product of their z-scored relative abundance level, representing the species-wise interactions for each subject.

After calculating the species-wise interactions, the modular structure obtained using the classic Louvain method³⁹ was employed to calculate subject-wise interaction strengths of each module by averaging all edges in a given module.

Resilience of microbial network

To determine the resilience of the network structure, simulated random and targeted attacks were conducted based on the binarized adjacent matrix. Random attack was performed by randomly removing one node and all its connections, while targeted attack was conducted by removing nodes and all its connections in a descending order based on the nodal degree. The influence of the simulated attacks to global efficiency of the network structure was evaluated after each removal. The global efficiency of node can be calculated as

2

where is the number of nodes within the network structure , and is the minimum path length between node and all the other nodes in the network. The global efficiency of the whole network was obtained by averaging across all nodes.

Calculating co-occurrence and co-exclusion interaction pairs of species

By calculating the sums for co-occurrence pairs and the differences for co-exclusion pairs, we generated metrics to represent the dynamics of species interactions within each pair. Specifically, co-occurrence pairs referred to species pairs that exhibited positive correlations, meaning their relative abundances tended to increase or decrease together. For each pair, we summed the relative abundances of the two species to create a single metric that reflects their combined abundance. In contrast, co-exclusion pairs were the species pairs that exhibited negative correlations, meaning when the relative abundance of one species increased, the other decreased. For these pairs, we calculated the absolute difference between the relative abundances of the two species. This difference quantifies the degree of imbalance or exclusivity in their abundances. These metrics were then used to explore potential associations between the interaction patterns and measures of cognition or temperament.

Extracting subject-wise trajectory features

Remarkable variations of microbiota trajectories not only among species but also subjects were observed. To characterize the patterns of different trajectories among subjects, we developed a clustering algorithm to group subjects into different clusters based on subject-wise trajectory patterns. Specifically, by examining the trajectories among subjects for a given microbe, one could visually separate subjects into two different groups; one group had lower relative abundance close to 0 while the other had higher relative abundance values. In addition, many subject-wise relative abundance trajectories may switch between low and high relative abundance groups during infancy. Moreover, the traditional zero-inflated model does not apply here because low-values were not exactly 0, and a thresholding approach could be biased by the choice of a threshold. Therefore, we developed the following two-step approach to model the data and perform clustering method on the trajectories. Specifically, a Gaussian mixture model was first used to model the distributions of the low level and high level abundance samples. The estimates of the Gaussian mixture model for each subject were extracted as features. Subsequently, a hierarchical clustering approach^54,55 using the complete-linkage method was used to cluster the subject-wise trajectories based on the features estimated from the Gaussian mixture model. The details of the above steps are provided below.

Gaussian mixture modeling

Let be the relative abundance for a given species of subject i at j-th visit, and is the age at the visit. We modeled the distribution of using the following Gaussian mixture model

3

4

5

6

where , are the regression coefficients for the corresponding model, is the random errors, and are the corresponding subject-wise random effects in the regression model and and indicate the average abundance level across multiple visits in the low-level and high-level groups for subject i, respectively. measures the probability of the samples of a subject belonging to the high-level group.

Hierarchical clustering

After obtaining the , and for each subject, we performed hierarchical clustering^54,55 using these three variables as features. Specifically, we first computed the Euclidean distance of the three features between subjects. The complete-linkage method of hierarchical clustering was then used to cluster the distance matrix and separate the subjects into different clusters. The average Silhouette distance was used to determine the optimal number of clusters.

Simulations

To evaluate if the number of visits of a given subject would affect the performance of the proposed Gaussian mixture model and the clustering approach, simulations were conducted. The number of subjects, number of visits per subject, and the age at each visit of the subjects were set to be the same as the study population. We randomly divided all subjects into 4 groups (Fig. S12, upper panel). In addition, the parameters for Eqs. 3–6 were set to . In addition, the random effects in the 4 groups were derived from the following distributions:

• Group 1:

• Group 2:

• Group 3:

• Group 4:

The settings above reflected that Groups 1 and 2 had trajectories in the high and low abundance groups, but Group 1 had a lower probability of transition to the high abundance group while Group 2 had a higher probability than the population. Groups 3 and 4 had a probability of being in the high abundance group. However, the abundance level in the lower abundance group was higher in group 3 and lower in group 4, and the opposite for the high abundance trajectories. The scatter plots of the simulated datasets were shown in Fig. S12.

We generated 100 independent datasets using the model above and the scatter plots of one simulated dataset were shown in Fig. S12. The predicted clustering labels using subjects with at least k visits, where k ranged between 1 and 5. The Cramer’s V statistic between the predicted labels and true labels were calculated to measure the accuracy of the proposed method. The results were shown in the lower panel of Fig. S12.

Statistical analysis

All statistical analysis were performed using R 3.6.1. Linear mixed effect models were used in our study to determine the potential associations between features of gut microbiota and cognition/temperament. Specifically, a random intercept was introduced for each subject, thus modeling the correlation of the observations from the same subjects at different visits. Moreover, demographic information including feeding practice, maternal education, and mode of delivery were controlled as covariates when either IBQ-R or ECBQ was used the outcome measures (Table S1). The Bonferroni correction was used for multiplicity adjustment with a significance level at 0.05. The linear mixed effect model was

7

where is the covariates of interest, is the random intercept for each subject, and is the random error. In contrast, since sex was shown to be significantly different for the MSEL composite as well as multiple subdomain scores (Table S2a), we included sex as a covariate in the model when MSEL was used as the outcome measures. Therefore, the model was

8

To determine the potential associations between the diversity measures of the gut microbiota and cognition/temperament, the age effects of diversity measures as well as temperament traits were regressed out using generalized additive mixed models (GAMM) before the analysis. The residuals (Figure S3) were then used for the aforementioned linear mixed effect models to discern the potential associations. The GAMM model was where is a cubic spline function with k knots, is the random intercept for each subject, and is the random error. We ranged k from 3 to 20, and the optimal k was chosen based on smallest AIC. The residuals were calculated as .

Acknowledgements

This study was supported in part by NIH grants (U01MH110274 (Lin, Elison); R01MH116527 (Li); R01MH086633 (Zhu); and MH1044324-03S1 (Elison)); MH015755 (Howell trainee); and a grant (Lin) from Nestlé Product Technology Center-Nutrition, Société des Produits Nestlé S.A., Switzerland. We thank all parents for consenting their children to participate in the study, Monika Chen, Bernard Berger and Léa Siegwald for help with organizing the microbiota analyses.

Author contributions

Designed research: BH, HH, JE, WL. Conducted research: ZZ, YY, BH, HH, JE, WL. Analyzed data or performed statistical analysis: ZZ, YY, TL, WY, SC, HZ, WL. Wrote paper: ZZ, YY, TMS, TL, WY, SC, HZ, NS. Had primary responsibility for final content: WL.

Data availability

The cognitive and temperament data have been uploaded to the NIMH Data Archive (https://nda.nih.gov/) as a part of the Baby Connectome Project. The microbiota data are available at European Nucleotide Archive (https://www.ebi.ac.uk/ena/browser/home) with the accession number of PRJEB73364 and Secondary Accession ERP158163.

Declarations

Competing interests

WL is a consultant of and had received travel support from Nestlé SA, Switzerland. TS, and NS are employees of Société des Produits Nestlé SA, Switzerland at the Nestlé Product Technology Center and Nestlé Institute of Health Sciences, respectively. ZZ is currently an employee of Google Inc, Mountain View CA, USA although the work reported in this manuscript was done when ZZ was a graduate student at UNC.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Jena, A et al. Gut-brain axis in the early postnatal years of life: A developmental perspective. Front. Integr. Neurosci.; 2020; 14, 44.1:CAS:528:DC%2BB3cXisVCrurzJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32848651][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7419604]

2. Collins, SM; Surette, M; Bercik, P. The interplay between the intestinal microbiota and the brain. Nat. Rev. Microbiol.; 2012; 10, pp. 735-742.1:CAS:528:DC%2BC38Xhtlylt7%2FJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23000955]

3. Cryan, JF; Dinan, TG. Mind-altering microorganisms: the impact of the gut microbiota on brain and behaviour. Nat. Rev. Neurosci.; 2012; 13, pp. 701-712.1:CAS:528:DC%2BC38XhtlamurvN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22968153]

4. Wang, S et al. Targeting the gut microbiota to influence brain development and function in early life. Neurosci. Biobehav. Rev.; 2018; 95, pp. 191-201. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30195933]

5. Cryan, JF et al. The microbiota-gut-brain axis. Physiol Rev; 2019; 99, pp. 1877-2013.1:CAS:528:DC%2BB3cXitFWjsLs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31460832]

6. Hui, Y et al. The effect of early probiotic exposure on the preterm infant gut microbiome development. Gut Microbes; 2021; 13, 1951113. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34264803][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8284123]

7. The infant gut microbiome as a microbial organ influencing host well-being. Italian J. Pediatr. https://doi.org/10.1186/s13052-020-0781-0.

8. De Angelis, M et al. Fecal microbiota and metabolome of children with autism and pervasive developmental disorder not otherwise specified. PLoS ONE; 2013; 8, 2013PLoSO..876993D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24130822][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3793965]e76993.

9. Inoue, R et al. A preliminary investigation on the relationship between gut microbiota and gene expressions in peripheral mononuclear cells of infants with autism spectrum disorders. Biosci Biotechnol Biochem; 2016; 80, pp. 2450-2458.1:CAS:528:DC%2BC28XhsVKjur%2FJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27581276]

10. Finegold, SM et al. Pyrosequencing study of fecal microflora of autistic and control children. Anaerobe; 2010; 16, pp. 444-453.1:CAS:528:DC%2BC3cXhtVymt73F [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20603222]

11. Kang, D-W et al. Reduced incidence of Prevotella and other fermenters in intestinal microflora of autistic children. PLoS ONE; 2013; 8, 2013PLoSO..868322K1:CAS:528:DC%2BC3sXhtFKkt7%2FF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23844187][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3700858]e68322.

12. Johnson, K. V.-A. Gut microbiome composition and diversity are related to human personality traits. Hum Microb J15, (2020).

13. Sherwin, E., Bordenstein, S. R., Quinn, J. L., Dinan, T. G. & Cryan, J. F. Microbiota and the social brain. Science366, (2019).

14. Brothers, L. Brain mechanisms of social cognition. J Psychopharmacol; 1996; 10, pp. 2-8.1:STN:280:DC%2BC383gtVClsQ%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22302721]

15. Roswall, J et al. Developmental trajectory of the healthy human gut microbiota during the first 5 years of life. Cell Host Microbe; 2021; 29, pp. 765-776.e3.1:CAS:528:DC%2BB3MXotVWhurY%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33794185]

16. Borre, YE et al. Microbiota and neurodevelopmental windows: Implications for brain disorders. Trends Mol Med; 2014; 20, pp. 509-518. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24956966]

17. Silbereis, JC; Pochareddy, S; Zhu, Y; Li, M; Sestan, N. The cellular and molecular landscapes of the developing human central nervous system. Neuron; 2016; 89, pp. 248-268.1:CAS:528:DC%2BC28XhtlKrtL8%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26796689][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4959909]

18. Thompson, RA; Nelson, CA. Developmental science and the media early brain development. Am Psychol; 2001; 56, pp. 5-15.1:STN:280:DC%2BD3M3hsFWktQ%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11242988]

19. Putnam, SP; Gartstein, MA; Rothbart, MK. Measurement of fine-grained aspects of toddler temperament: The early childhood behavior questionnaire. Infant Behav Dev; 2006; 29, pp. 386-401. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17138293][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4334385]

20. Temperament, personality and developmental psychopathology: A review based on the conceptual dimensions underlying childhood traits. Child Psychiatry Human Dev. https://doi.org/10.1007/s10578-009-0171-8.

21. Compas, BE; Connor-Smith, J; Jaser, SS. Temperament, stress reactivity, and coping:implications for depression in childhood and adolescence. J Clin Child Adolesc Psychol; 2004; 33, pp. 21-31. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15028538]

22. De Pauw, SSW; Mervielde, I. The role of temperament and personality in problem behaviors of children with ADHD. J Abnorm Child Psychol; 2011; 39, pp. 277-291. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20862537]

23. Izard, CE; Schultz, D; Fine, SE; Youngstrom, E; Ackerman, BP. Temperament, cognitive ability, emotional knowledge, and adaptive social behavior. Imagin. Cogn. Pers.; 1999; 19, pp. 305-330.

24. Contextual risk factors impacting the colonization and development of the intestinal microbiota: Implications for child…. http://ouci.dntb.gov.ua/en/works/9QARGYJ9/.

25. Aatsinki, A-K et al. Gut microbiota composition is associated with temperament traits in infants. Brain Behav. Immun.; 2019; 80, pp. 849-858. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31132457]

26. Carlson, AL et al. Infant gut microbiome associated with cognitive development. Biol Psychiatry; 2018; 83, pp. 148-159. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28793975]

27. Christian, LM et al. Gut microbiome composition is associated with temperament during early childhood. Brain Behav. Immun.; 2015; 45, pp. 118-127. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25449582]

28. Loughman, A et al. Gut microbiota composition during infancy and subsequent behavioural outcomes. EBioMedicine; 2020; 52, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32062351][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7016366]102640.

29. Rammstedt, B; Lechner, CM; Danner, D. Relationships between personality and cognitive ability: A facet-level analysis. J. Intell.; 2018; 6, 28. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31162455][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6480763]

30. Fox, M et al. Development of the infant gut microbiome predicts temperament across the first year of life. Dev Psychopathol; 2021; [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36914287][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8916713][DOI: https://dx.doi.org/10.1017/S0954579421000456]

31. Stewart, CJ et al. Temporal development of the gut microbiome in early childhood from the TEDDY study. Nature; 2018; 562, pp. 583-588.2018Natur.562.583S1:CAS:528:DC%2BC1cXitVSitLnN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30356187][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6415775]

32. Mullen, E. Mullen scales of early learning. American Guidance Service (1995).

33. Gartstein, MA; Rothbart, MK. Studying infant temperament via the revised infant behavior questionnaire. Infant. Behav. Dev.; 2003; 26, pp. 64-86.

34. Sordillo, JE et al. Association of the infant gut microbiome with early childhood neurodevelopmental outcomes: an ancillary study to the VDAART randomized clinical trial. JAMA Netw Open; 2019; 2, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30901046][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6583279]e190905.

35. Zelezniak, A et al. Metabolic dependencies drive species co-occurrence in diverse microbial communities. Proc. Natl. Acad. Sci. USA; 2015; 112, pp. 6449-6454.2015PNAS.112.6449Z1:CAS:528:DC%2BC2MXns1Oks78%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25941371][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4443341]

36. Hall, CV et al. Co-existence of network architectures supporting the human gut microbiome. iScience; 2019; 22, pp. 380-391.2019iSci..22.380H [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31812808][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6911941]

37. Friedman, J; Alm, EJ. Inferring correlation networks from genomic survey data. PLoS Comput. Biol.; 2012; 8, 2012PLSCB..8E2687F1:CAS:528:DC%2BC38XhsVaktr3M [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23028285][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3447976]e1002687.

38. Weiss, S et al. Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. ISME J.; 2016; 10, pp. 1669-1681.1:CAS:528:DC%2BC28XhtVahu7fF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26905627][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4918442]

39. Blondel, VD; Guillaume, J-L; Lambiotte, R; Lefebvre, E. Fast unfolding of communities in large networks. J. Stat. Mech. Theory Exp.; 2008; P10008, pp. 1-12.

40. Faskowitz, J; Esfahlani, FZ; Jo, Y; Sporns, O; Betzel, RF. Edge-centric functional network representations of human cerebral cortex reveal overlapping system-level architecture. Nat. Neurosci.; 2020; 23, pp. 1644-1654.1:CAS:528:DC%2BB3cXitFWju7nM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33077948]

41. Strati, F et al. New evidences on the altered gut microbiota in autism spectrum disorders. Microbiome; 2017; 5, 24. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28222761][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5320696]

42. Ren, T et al. Gut microbiota altered in mild cognitive impairment compared with normal cognition in sporadic Parkinson’s disease. Front. Neurol.; 2020; 11, 137. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32161568][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7052381]

43. Baxter, NT et al. Dynamics of human gut microbiota and short-chain fatty acids in response to dietary interventions with three fermentable fibers. MBio; 2019; 10, pp. e02566-e2618.1:CAS:528:DC%2BC1MXht12htL7I [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30696735][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6355990]

44. Ze, X; Duncan, SH; Louis, P; Flint, HJ. Ruminococcus bromii is a keystone species for the degradation of resistant starch in the human colon. ISME J; 2012; 6, pp. 1535-1543.1:CAS:528:DC%2BC38XhtVOjtL3L [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22343308][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3400402]

45. Eastwood, J; Walton, G; Van Hemert, S; Williams, C; Lamport, D. The effect of probiotics on cognitive function across the human lifespan: A systematic review. Neurosci. Biobehav. Rev.; 2021; 128, pp. 311-327. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34171323]

46. Firmansyah, A et al. Improved growth of toddlers fed a milk containing synbiotics. Asia Pac. J. Clin. Nutr.; 2011; 20, pp. 69-76.1:CAS:528:DC%2BC3MXltlyrsbc%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21393113]

47. Sudo, N et al. Postnatal microbial colonization programs the hypothalamic–pituitary–adrenal system for stress response in mice. J. Physiol.; 2004; 558, pp. 263-275.1:CAS:528:DC%2BD2cXmtFWhtro%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15133062][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1664925]

48. Diaz-Heijtz, R et al. Normal gut microbiota modulates brain development and behavior. Proc. Natl. Acad. Sci. USA; 2011; 108, pp. 3047-3052.2011PNAS.108.3047H [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21282636]

49. Desbonnet, L; Clarke, G; Shanahan, F; Dinan, TG; Cryan, JF. Microbiota is essential for social development in the mouse. Mol Psychiatry; 2014; 19, pp. 146-148.1:STN:280:DC%2BC3snlt1aksg%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23689536]

50. Su, PL; Rogers, SJ; Estes, A; Yoder, P. The role of early social motivation in explaining variability in functional language in toddlers with autism spectrum disorder. Autism; 2021; 25, pp. 244-257. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32921137]

51. Mulvey, N; Jenkins, L. Language skills as predictors of social skills and behaviors in preschool children. Contemp. Sch. Psychol.; 2021; 25, pp. 503-514.

52. Microbial community profiling of human saliva using shotgun metagenomic sequencing - PubMed. https://pubmed.ncbi.nlm.nih.gov/24846174/.

53. Watts, SC; Ritchie, SC; Inouye, M; Holt, KE. FastSpar: rapid and scalable correlation estimation for compositional data. Bioinformatics; 2019; 35, pp. 1064-1066.1:CAS:528:DC%2BC1MXitVWgsrjP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30169561]

54. Data mining and knowledge discovery handbook. (Springer, New York, 2005). https://doi.org/10.1007/b107408.

55. Nielsen, F. Introduction to HPC with MPI for Data Science; 2016; Cham, Springer International Publishing: [DOI: https://dx.doi.org/10.1007/978-3-319-21903-5]