Microbiome function and neurodevelopment in Black infants: vitamin B₁₂ emerges as a key factor

ABSTRACT

The early life gut microbiome affects the developing brain, and therefore may serve as a target to support neurodevelopment of children living in stressful and under-resourced environments, such as Black youth living on the South Side of Chicago, for whom we observe racial disparities in health. Microbiome compositions/functions key to multiple neurodevelopmental facets have not been studied in Black children, a vulnerable population due to racial disparities in health; thus, a subsample of Black infants living in urban, low-income neighborhoods whose mothers participated in a prenatal nutrition study were recruited for testing associations between composition and function of the gut microbiome (16S rRNA gene sequencing, shotgun metagenomics, and targeted metabolomics of fecal samples) and neurodevelopment (developmental testing, maternal report of temperament, and observed stress regulation). Two microbiome community types, defined by high Lachnospiraceae or Enterobacteriaceae abundance, were discovered in this cohort from 16S rRNA gene sequencing analysis; the Enterobacteriaceae-dominant community type was significantly negatively associated with cognition and language scores, specifically in male children. Vitamin B₁₂ biosynthesis emerged as a key microbiome function from shotgun metagenomics sequencing analysis, showing positive associations with all measured developmental skills (i.e., cognition, language, motor, surgency, effortful control, and observed stress regulation). Blautia spp. also were identified as substantial contributors of important microbiome functions, including vitamin B₁₂ biosynthesis and related vitamin B₁₂-dependent microbiome functions, anti-inflammatory microbial surface antigens, competitive mechanisms against pathobionts, and production of antioxidants. The results are promising with respect to the potential for exploring therapeutic candidates, such as vitamin B₁₂ nutritional or Blautia spp. probiotic supplementation, to support the neurodevelopment of infants at risk for experiencing racial disparities in health.

KEYWORDS:

• Neurodevelopment
• Microbiome
• Vitamin B₁₂
• Black infant
• Health disparities

Introduction

Low resources and high stress negatively impact cognitive, emotional, and behavioral development in children.^1,2 There is a pressing need to identify modifiable targets that can reduce these disparities in neurodevelopmental outcomes, of which the gut microbiome (i.e., microorganisms that inhabit the gastrointestinal tract) shows promise. The gut microbiome begins as a small species consortium that diversifies as the host matures; vertical transmission of microbes from the mother after birth is a significant source of the earliest colonizers in the infant gut.³ The gut microbiome plays vital roles in its host physiology that affect the developing brain by contributing to digestion followed by producing neuroactive fermentative by-products, synthesizing micronutrients, and modulating the immune system.^4,5 Both compositional and functional deviations of the microbiome (i.e., dysbiosis) have been directly associated with deficits in attention, and social and emotional regulation.^6,7 Living in an under-resourced environment is associated with gut microbiome dysbiosis,^8,9 and a dysbiotic microbiome can be passed on from family to child,¹⁰ possibly contributing to the persistence of health disparities. Microbiome-based therapeutics could contribute to reducing the severity of or removing such health disparities. However, it is not understood which gut microbiome colonizers nor which microbiome functions in early life are critical to neurodevelopment for children developing in stressful and under-resourced social environments.

This study is the first to thoroughly examine the microbiome at both a compositional and functional level in the context of multiple aspects of neurodevelopment simultaneously for a cohort of Black children whose families live in under-resourced environments. We leveraged data from an existing cohort of infants whose mothers had participated in a randomized controlled trial of omega-3 fatty acid supplementation during pregnancy (NAPS: NCT02647723). The Nutrition and Pregnancy Study (NAPS) aimed to use nutrition to improve prenatal stress regulation in Black women exposed to high levels of unpredictable stressors (e.g., community violence), resource-related stressors (e.g., trouble finding affordable housing), and discrimination stress as a means of improving maternal health, and birth and neurodevelopmental outcomes. The conceptual framework for the NAPS Study was grounded in the Developmental Origins of Health and Disease and Prenatal Programming models of health.¹¹ In this context specifically, we posited that a primary cause of racial disparities in maternal and child health among Black Americans begins with high exposure to psychosocial stress during pregnancy. Growing evidence indicates that DHA supplementation improves maternal stress reactivity during pregnancy and protects neurodevelopment of the offspring, especially in the context of high levels of stress exposure.^12,13 In the NAPS study, we aimed to translate such findings into a nutrition-based preventive intervention program to reduce health disparities for Black women and children. Women were eligible for the trial if they identified as Black, were insured through Medicaid and received prenatal care at the University of Chicago or Federally Qualified Health Center affiliates. Medicaid is a joint state and federal government health insurance program in the United States for individuals and families with low-incomes.¹⁴ Infants were seen at 4- and 9-months for assessments of multiple functional domains, including verbal, motor and cognitive development through standardized testing, temperament from maternal report, and observed stress regulation determined by a laboratory-based probe (i.e., cortisol reactivity after a stressor). Infant fecal samples were collected from a subsample of NAPS participants (n = 28; 1 sample per participant) that were age-eligible at the time of study approval and subjected to 16S rRNA gene and shotgun metagenomics sequencing to thoroughly interrogate the composition and potential functions of the gut microbiome. Linear regression models adjusted for gestational age at birth, birthweight, and age at assessment revealed significant associations between richness, Shannon diversity, and overall β-diversity and multiple indices of neurodevelopment. Further, the significant associations were moderated by sex such that associations between microbiome composition and neurodevelopment was observed primarily among male infants. The gut microbiome function of vitamin B₁₂ biosynthesis was significantly associated with multiple indices in the developmental skills domain (i.e., cognition, language, motor, surgency, effortful control, cortisol reactivity after a stressor) in male children as were several vitamin B₁₂-dependent metabolic pathways. Together, these results indicate a cascading effect of vitamin B₁₂ depletion that resulted in dysregulated microbiome metabolic output and pathogenicity.

Results

Cohort characteristics

The goal of this study was to examine the associations between the infant gut microbiome (n = 28; 1 sample per participant) and neurodevelopment. The 28 participants in the present study are children whose mothers were enrolled in a randomized controlled trial (RCT) of prenatal docosahexaenoic acid (DHA) supplementation aimed at improving maternal stress regulation during pregnancy (NAPS: NCT02647723). The full sample for the RCT was 168 women; the data are currently being analyzed, and thus results from the RCT are not available; of the 28 children in the present study, 20 were born to mothers who were assigned to the active supplement arm and 8 to mothers assigned to the placebo arm. For the purposes of transparency and generalizability of the results from the present study, we compared maternal DHA blood levels at 36 weeks of gestation between the 20 children whose mothers received DHA supplementation and the 8 children whose mothers received placebo; the comparison was not statistically significant (Table S1).

We used the Great Cities Community Area Hardship Index,^15,16 which averages scores for six variables including unemployment, education, per capita income level, poverty, crowded housing, and dependency. The range for the entire Chicago Community is 9.4 to 76.5, wherein higher scores indicate worse economic conditions. For the 28 participants in the present study, hardship scores ranged from 36.5 to 76.5, with a median index of 52.9; neighborhoods with a hardship index score of 52.9 are characterized as having an unemployment rate of 19% and a per capita income of 20,782 USD.

Additional demographic characteristics of the participants and descriptive statistics for the neurodevelopmental assessments are presented in Table 1. The majority of participants were vaginally delivered, term infants of average birthweight with an even distribution of males (n = 14) and females (n = 14). The average age at assessment for the present study was 7 months (range = 3–16 months). Delivery mode was not significantly associated with any of the neurodevelopmental metrics, but birthweight was significantly (p < .05) negatively associated with BSID-III motor composite and ICQ fussy/difficult scores (Table S2). Further, compared to female infants, male infants had significantly lower BSID-III cognitive and language composite scores (Table S2). Since birthweight is resultant of the in utero environment and thus not influenced by the child’s microbiome, this variable was adjusted for in all statistical models, along with gestational age at birth plus age at assessment visit which knowingly impact the results of neurodevelopmental assessments. However, the interaction effects of sex with the microbiome were of interest to determine, and further investigation of how the study covariates differed by sex was performed.

Table 1. Child cohort demographics and neurodevelopmental outcome metrics.

Demographics	Mean ± standard deviation
Gestational age at birth (weeks)	39.17 ± 1.11
Birthweight (kg)	3.088 ± 0.381
Age at evaluation visit (months old)	7.59 ± 3.65
Sex, % male	50%
Delivery mode, % vaginal delivery	81%
Feeding type, % formula (4 months old)	95%
Feeding type, % formula (9 months old)	95%
Solid foods, % in last 7 days (4 months old)	37%
Solid foods, % in last 7 days (9 months old)	58%
Amount of night sleep, hours (4 months old)	7.88 ± 1.82
Amount of night sleep, hours (9 months old)	8.53 ± 1.58
Amount of day sleep, hours (4 months old)	5.45 ± 3.01
Amount of day sleep, hours (9 months old)	4.58 ± 2.55
Number of night wakenings (4 months old)	1.65 ± 0.81
Number of night wakenings (9 months old)	1.00 ± 0.94
Time to put child to sleep, minutes (4 months old)	21.85 ± 13.43
Time to put child to sleep, minutes (9 months old)	27.58 ± 39.39
Bayley Scales of Infant Development III	Mean ± standard deviation
Cognitive composite score	101.90 ± 19.14
Language composite score	88.00 ± 13.02
Expressive communication scaled score	7.20 ± 2.10
Receptive communication scaled score	8.48 ± 3.34
Motor composite score	95.65 ± 15.19
Fine motor scaled score	9.125 ± 3.88
Gross motor scaled score	10.36 ± 3.46
Infant Behavior Questionnaire Revised	Mean ± standard deviation
Effortful control score	5.51 ± 1.01
Surgency score	5.23 ± 1.12
Negative affect score	3.82 ± 1.27
Infant Characteristics Questionnaire	Mean ± standard deviation
Dull score	7.46 ± 3.58
Fussy/Difficult score	14.67 ± 5.16
Unadaptable score	10.75 ± 4.96
Unpredictable score	5.75 ± 3.33
Cortisol response to Face-to-Face Still-Face paradigm	Mean ± standard deviation
Cortisol reactivity (AUCI)	−2.15 ± 9.32

Approximately 70% of study participants completed maternal reports of child feeding and sleep practices at 4- and 9-months of age. From this and earlier data, it was determined that delivery mode, birthweight, age at neurodevelopmental assessment visit, and all sleep metrics (amount of day sleep, amount of night sleep, number of nighttime wakenings, and time to put child to sleep) did not significantly differ by sex (Table S3). 95% of the infants were fed formula (Table 1). At the age 4-month assessment, male children were significantly (p < .05) more likely to have received solid foods in the last 7 days than female children (Table S3). This finding was no longer statistically significant at the age 9-month assessment (Table S3). Statistical analysis of the entire NAPS cohort (i.e., with and without microbiome data) to corroborate these findings (94 infants with data on feeding at 4-months and 80 infants with data on feeding at 9-months) again found that the only statistically significant difference (p < .05) between male and female children in the entire NAPS cohort was that male children were significantly more likely to receive solid foods in the past 7 days at 4 months of age (Table S3), thus verifying our findings in the studied subsample of participants with microbiome data.

Assessment of developmental skills, temperament, and observed stress regulation reveals two neurodevelopmental domains

Cluster analysis of the neurodevelopmental outcome metrics revealed that these assessments could be grouped into two neurodevelopmental domains (Figure 1). The first and largest domain was comprised of developmental skills in cognition, language and motor, temperamental traits reflecting engagement, and observed stress regulation (i.e., cortisol reactivity to a stressor), which we refer to as the ‘Developmental skills’ domain. The second domain encompassed temperamental traits reflecting difficulty with resilience and/or adaptation that we refer to as the ‘Difficult temperament’ domain.

Figure 1. Partial Pearson correlation of neurodevelopmental outcome metrics for a cohort of Black children living in low socioeconomic and high stress environments (average age 7.6 months old; n = 28), after adjustment for gestational age at birth, birthweight, age at assessment visit and sex. Metrics are ordered by hierarchical clustering via the ward linkage method, and optimal clusters are enclosed in rectangles as chosen by the permutation around medoids (i.e., k-medoids) clustering solution with the largest average silhouette width. Significant correlations are indicated as follows: ***p < .001; **p < .01; *p < .05.

Microbiome structure was associated with developmental skills in Black infants

Two microbiome community types were natively found in the NAPS infants, and these community types were specifically associated with the developmental skills domain in males. Fecal microbiome composition was first examined by shotgun metagenomics to understand the distribution of microbial kingdoms across individuals. The microbiome was comprised of 100% bacteria for all but two children for whom the microbiome contained small abundances of Candida. Thus, community typing analysis focused on bacteria, and the probabilistic method of Dirichlet Multinomial Mixtures was applied to the 16S rRNA gene sequencing data.

One community type was found at the taxonomic levels of phylum and class, indicating cohort conformity at the higher taxonomic levels. Two community types were found at the taxonomic levels of family and genus (Figure 2a). Community type 1 was differentiated by a high abundance of the bacterial family Lachnospiraceae, and the bacterial genera Blautia, Bifidobacterium and Bacteroides. Community type 2 was differentiated by a high abundance of the bacterial family Enterobacteriaceae, which were unclassified at the genus level. The shotgun metagenomics data was next examined to determine which species of these microbial taxa were elevated by abundance (>1% mean difference) and presence (>10% more prevalent) in the specific community types. The top species for each community type 1 genus were: Blautia hansenii, Blautia coccoides, and Blautia wexlerae; Bifidobacterium pseudocatenulatum and Bifidobacterium dentium; and Bacteroides vulgatus, Bacteroides uniformis and Bacteroides thetaiotaomicron (Table S4). The top species of Enterobacteriaceae in community type 2 were: Klebsiella variicola, Klebsiella pneumoniae, and Klebsiella quasipneumoniae (Table S4). Finally, Gene-Set Enrichment Analysis (GSEA) of the Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologies (KOs) was undertaken using the t-scores from linear models to uncover significant differences (p < .05; false-positive rate < 1%) in microbiome functional pathways between community types (Table S5). The top 5 metabolic pathways differentiating the community types were: Biosynthesis of amino acids, Biosynthesis of cofactors, Biosynthesis of various other secondary metabolites, Pyrimidine metabolism, and 2-Oxocarboxylic acid metabolism (Figure 2a).

Figure 2. Overall α-diversity and β-diversity of the fecal microbiome (n = 28) in relation to the neurodevelopmental outcome metrics. (a) Two community types were determined to be present in the child (average age 7.6 months old) fecal microbiome dataset through the probabilistic method of dirichlet multinomial mixtures. These community types could be described at the taxonomic levels of family and genus, and top microbial genera driving the differential composition of these community types are indicated. Statistically significant Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways differentiating the two community types from gene-set enrichment analysis of the t-scores produced from linear modeling are also displayed by the absolute value of their normalized enrichment scores (NES) and relative significance. ***p < .0001; **p < .001; *p < 0.01. (b) Overall α-diversity (i.e., richness) and β-diversity (i.e., community type) were significantly associated with the Bayley Scales of Infant Development III cognitive and language composite scores in a sex-dependent manner. Significance (p < .05) was determined through multiple regression after adjusting for gestational age at birth, birthweight, and age at assessment visit. Interaction plots of the standardized variables after modeling demonstrate the sex-dependent relationship of richness and these neurodevelopmental outcomes. Estimated marginal means (emmean) with 95% confidence intervals from the standardized variables after modeling demonstrate the sex-dependent relationship of community type and these neurodevelopmental outcomes.

We next tested the relationship between the microbiome community types and the two neurodevelopmental domains. Richness, Shannon diversity and community type 1 were significantly (p < .05) associated with the BSID-III language composite score, specifically receptive communication (Table S6). Additionally, both richness and community type 1 were significantly associated with the BSID-III cognitive composite score (Table S6). These effects were significantly moderated (p < .05) by sex: among males, an Enterobacteriaceae-dominated microbiome was associated with suboptimal neurodevelopment, whereas no association between an Enterobacteriaceae-dominated community type and neurodevelopment was found for females (Figure 2b). Importantly, no statistically significant differences were found for fecal microbiome α-diversity or community type between males and females (Table S6). Thus, despite equal representation of an Enterobacteriaceae-dominated community type for males and females, dysbiosis conferred risk to brain development and function only among male infants. This result warrants the inclusion of sex as a moderator for regression models, rather than segregating statistical analysis by sex.

Biosynthesis of vitamin B₁₂ is the key metabolic pathway driving associations between the developmental skills domain and the microbiome in males

Investigation of the significant metabolic pathways associated with each neurodevelopmental assessment revealed that biosynthesis of the cofactor cobalamin (vitamin B₁₂) was the most strongly associated pathway with the developmental skills domain in a sex-dependent manner. Several microbial pathways containing vitamin B₁₂-dependent enzymes were also significantly associated, suggesting a cascading effect of vitamin B₁₂ depletion on microbiome succession and functional output with consequences for the developing male host. To interrogate which pathways were specifically associated with the neurodevelopmental metrics, Gene-Set Enrichment Analysis (GSEA) on the t-scores from multiple regression models of the Kyoto Encyclopedia of Genes and Genomes (KEGG) orthologies (KOs) that adjusted for gestational age at birth, birthweight, age at assessment visit, and considering sex as a moderator were conducted for each neurodevelopmental assessment. The KEGG modules and pathways were then tallied to determine how many facets of each of the two neurodevelopmental domains were significantly (p < 0.05; false-positive rate < 1%) associated.

Vitamin B₁₂ biosynthesis was the only microbiome function significantly associated with all neurodevelopmental metrics in the developmental skills domain (i.e., cognition, language, motor, surgency, effortful control, cortisol reactivity; Table S7), and was the most significantly positively associated KEGG module by t-score for the BSID-III cognitive composite score, receptive communication scaled score, and motor composite score, including both fine and gross motor skills, in addition to being the most significantly negatively associated KEGG module by t-score for cortisol reactivity (Table 2). Whereas the associations of vitamin B₁₂ biosynthesis with motor skills and cortisol reactivity were sex-independent, the impacts on the other metrics in the developmental skills domain were exclusive to males (Figure 3; Table S8). Vitamin B₁₂ biosynthesis capabilities were found in many Blautia spp., in addition to several other Lachnospiraceae spp. associated with community type 1 (Figure 2a; Table 3). A few Bacteroides spp. were also able to synthesize vitamin B₁₂, including the species most strongly associated with community type 1, Ba. vulgatus (Figure 2a; Table 3).

Table 2. Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways and modules of the fecal microbiome that were most significantly positively and negatively associated with each neurodevelopmental outcome metric.

Developmental skills domain neuro-developmental metric	KEGG pathway and [module]	Normalized enrichment score (p value)
BSID-III cognitive composite score	Secondary bile acid biosynthesis ↑ [Cobalamin biosynthesis] ↑ Biosynthesis of amino acids ↓ [Deoxyribonucleotide biosynthesis] ↓	2.4 (p = 7.5 × 10^−6) 2.0 (p = 0.00023) −1.8 (p = 9.5 × 10^−6) −1.7 (p = 3.4 × 10^−6)
BSID-III language composite score	Lipoarabinomannan (LAM) biosynthesis ↑ [Methylaspartate cycle] ↑ Biofilm formation – Escherichia coli ↓ [None significant] ↓	2.3 (p = 9.7 × 10^−7) 2.1 (p = 0.00047) −2.1 (p = 3.1 × 10^−6) N/A
BSID-III expressive communication scaled score	Lipoarabinomannan (LAM) biosynthesis ↑ [Tryptophan metabolism] ↑ Flagellar assembly ↓ [Phenylacetate degradation] ↓	2.7 (p = 1.7 × 10^−6) 1.8 (p = 0.0070) −2.1 (p = 1.0 × 10^−7) −1.7 (p = 0.0041)
BSID-III receptive communication scaled score	Lipoarabinomannan (LAM) biosynthesis ↑ [Cobalamin biosynthesis] ↑ Biofilm formation – Escherichia coli ↓ [Coenzyme A biosynthesis] ↓	2.3 (p = 1.3 × 10^−6) 2.1 (p = 0.00030) −2.1 (p = 1.9 × 10^−5) −1.8 (p = 0.0028)
BSID-III motor composite score	Porphyrin metabolism ↑ [Cobalamin biosynthesis] ↑ Pathogenic Escherichia coli infection ↓ [Phenylacetate degradation] ↓	2.1 (p = 4.3 × 10^−5) 2.1 (p = 0.00035) −1.9 (p = 0.00045) −1.7 (p = 0.0068)
BSID-III fine motor scaled score	O-Antigen nucleotide sugar biosynthesis ↑ [Cobalamin biosynthesis] ↑ Biosynthesis of amino acids ↓ [Leucine biosynthesis] ↓	2.0 (p = 0.00023) 2.0 (p = 0.00058) −2.1 (p = 1.5 × 10^−8) −1.7 (p = 0.0013)
BSID-III gross motor scaled score	Teichoic acid biosynthesis ↑ [Cobalamin biosynthesis] ↑ Biofilm formation – Escherichia coli ↓ [None significant] ↓	2.1 (p = 0.00036) 2.2 (p = 0.0014) −1.9 (p = 2.0 × 10^−5) N/A
IBQ-R effortful control score	Lipoarabinomannan (LAM) biosynthesis ↑ [Methylaspartate cycle] ↑ Biosynthesis of amino acids ↓ [De novo pyrimidine biosynthesis] ↓	2.0 (p = 0.00061) 1.9 (p = 0.0014) −2.3 (p = 1 × 10^−10) −2.1 (p = 2.7 × 10^−5)
IBQ-R surgency score	Lipoarabinomannan (LAM) biosynthesis ↑ [Methanogenesis] ↑ Biosynthesis of amino acids ↓ [De novo pyrimidine biosynthesis] ↓	2.0 (p = 8.9 × 10^−5) 1.8 (p = 0.0016) −2.1 (p = 3.1 × 10^−8) −2.0 (p = 0.00046)
Cortisol reactivity (AUCI)	Phenylalanine metabolism ↑ [Ethylmalonyl pathway] ↑ None significant ↓ [Cobalamin biosynthesis] ↓	2.0 (p = 0.00013) 1.7 (p = 0.0013) N/A −1.9 (p = 0.0050)
Difficult temperament domain neuro-developmental metric	KEGG pathway and [module]	Normalized enrichment score (p value)
IBQ-R negative affect score	Lipopolysaccharide biosynthesis ↑ [KDO2-lipid A biosynthesis] ↑ Lipoarabinomannan (LAM) biosynthesis ↓ [Purine degradation] ↓	2.2 (p = 2.2 × 10^−5) 2.2 (p = 1.1 × 10^−5) −2.3 (p = 8.9 × 10^−6) −2.0 (p = 0.0012)
ICQ dull score	Biosynthesis of amino acids ↑ [De novo pyrimidine biosynthesis] ↑ Steroid degradation ↓ [Dissimilatory sulfate reduction] ↓	2.0 (p = 1.6 × 10^−7) 2.0 (p = 0.00033) −1.7 (p = 0.00028) −1.8 (p = 0.00036)
ICQ fussy/difficult score	Histidine metabolism ↑ [Coenzyme A biosynthesis] ↑ Lipoarabinomannan (LAM) biosynthesis ↓ [Trans-cinnamate degradation] ↓	2.3 (p = 0.00015) 1.9 (p = 0.0067) −1.9 (p = 2.8 × 10^−5) −1.7 (p = 0.0072)
ICQ unadaptable score	Bacterial chemotaxis ↑ [Lysine biosynthesis] ↑ Benzoate degradation ↓ [Purine degradation] ↓	1.8 (p = 0.0091) 1.9 (p = 0.0012) −2.2 (p = 3.3 × 10^−6) −2.0 (p = 9.9 × 10^−5)
ICQ unpredictable score	Histidine metabolism ↑ [Histidine biosynthesis] ↑ Flagellar assembly ↓ [Denitrification] ↓	2.3 (p = 7.3 × 10^−5) 2.6 (p = 1.1 × 10^−5) −2.2 (p = 5.8 × 10^−6) −1.8 (p = 0.00071)

Figure 3. Differential Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways or modules of the fecal microbiome associated with child (average age 7.6 months old; n = 28) neurodevelopmental outcome metrics by sex after gene-set enrichment analysis of t-scores generated from multiple regression models that adjusted for gestational age at birth, birthweight, and age at assessment visit. Estimated marginal means (emmean) for a representative KEGG orthology from each of these pathways/modules are displayed for each sex, with magnitude indicated by the size of circle and direction indicated by the color (yellow = positive association; purple = negative association). Significance (p < .05; false-positive rate < 1%) is indicated from either the model regression term (red colored circle border) or sex interaction (green colored circle border).

Table 3. Microbial taxa contributing to the fecal microbiome Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways/modules significantly associated with the neurodevelopmental domains.

Vitamin B12 associated KEGG pathway/module	Contributing microbial taxa
Cobalamin biosynthesis	Bacteroides dorei Bacteroides vulgatus Blautia spp. Citrobacter freundii Clostridioides difficile Coprococcus catus Corynebacterium variabile Escherichia coli Eubacterium hallii Faecalibacterium prausnitzii Flavonifractor plautii Klebsiella spp. Lactococcus lactis Roseburia hominis Roseburia intestinalis Ruminococcus bicirculans Streptococcus parasanguinis Veillonella dispar Veillonella parvula
Deoxyribonucleotide biosynthesis/De novo pyrimidine biosynthesis (nrdJ exclusively)	Acidaminococcus intestini Bacteroides stercoris Blautia hansenii Corynebacterium falsenii Lactobacillus ruminis
Deoxyribonucleotide biosynthesis/De novo pyrimidine biosynthesis (nrdAB/nrdEF)	Bifidobacterium spp. Blautia spp. Citrobacter spp. Corynebacterium variabile Enterobacter cloacae complex Escherichia coli Haemophilus parainfluenzae Haemophilus sp. HMSC71H05 Klebsiella spp. Leuconostoc citreum Paeniclostridium sordellii Staphylococcus aureus Staphylococcus epidermidis
Ethylmalonyl pathway	Alistipes shahii Citrobacter freundii Escherichia coli Klebsiella spp. Odoribacter splanchnicus
Methylaspartate cycle	Blautia spp. Citrobacter freundii Escherichia coli Klebsiella spp.
Methanogenesis	Blautia coccoides Blautia producta Eubacterium callanderi Eubacterium limosum
Microbial antigen and pathogenicity associated KEGG pathway/module	Contributing microbial taxa
Lipoarabinomannan (LAM) biosynthesis	Bifidobacterium animalis Corynebacterium falsenii Corynebacterium variabile
Teichoic acid biosynthesis	Blautia coccoides Lactobacillus spp. Lactococcus lactis Lactococcus petauri Leuconostoc citreum Paeniclostridium sordellii Streptococcus spp.
Pathogenic E. coli infection	Escherichia coli
Bacterial chemotaxis	Escherichia coli Klebsiella oxytoca Klebsiella pneumoniae
Biofilm formation E. coli	Escherichia coli Klebsiella oxytoca Klebsiella pneumoniae
KDO2-lipid A biosynthesis	Acidaminococcus intestini Alistipes finegoldii Alistipes shahii Bacteroides spp. Escherichia coli Haemophilus parainfluenzae Klebsiella oxytoca Klebsiella pneumoniae Veillonella atypica
Flagellar assembly	Eubacterium eligens Eubacterium rectale Flavonifractor plautii Klebsiella pneumoniae Roseburia intestinalis
Secondary bile acid biosynthesis (7α-dehydroxylation)	Clostridium scindens Dorea sp. D27
Secondary bile acid biosynthesis (3α-dehydrogenation)	Eggerthella lenta Gordonibacter pamelaeae
dTDP-L-rhamnose biosynthesis	Alistipes shahii Anaerostipes hadrus Bacteroides spp. Bifidobacterium longum Blautia spp. Clostridium sp. CAG-299 Coprococcus catus Corynebacterium falsenii Corynebacterium variabile Eubacterium hallii Klebsiella michiganensis Klebsiella pneumoniae Lactobacillus ruminis Odoribacter splanchnicus
Valine, leucine and isoleucine degradation	Acidaminococcus intestini Anaerostipes hadrus Bacteroides spp. Bifidobacterium animalis Blautia spp. Candida tropicalis Coprococcus eutactus Corynebacterium falsenii Corynebacterium variabile Eggerthella lenta Escherichia coli Gordonibacter pamelaeae Haemophilus parainfluenzae Haemophilus sp. HMSC71H05 Klebsiella oxytoca Klebsiella pneumonia Lawsonella clevelandensis Roseburia spp. Rothia mucilaginosa Streptococcus parasanguinis
Purine, histidine and phenylalanine associated KEGG pathway/module	Contributing microbial taxa
Purine degradation	Blautia spp. Coprococcus catus Corynebacterium variabile Enterobacter cloacae complex Escherichia coli Klebsiella oxytoca Klebsiella pneumoniae
Histidine biosynthesis	Acidaminococcus intestini Alistipes finegoldii Anaerostipes hadrus Bacteroides spp. Bifidobacterium spp. Blautia coccoides Ruminococcus torques Citrobacter freundii Clostridioides difficile Coprococcus catus Corynebacterium falsenii Corynebacterium variabile Enterobacter cloacae complex Escherichia coli Eubacterium eligens Klebsiella spp. Lactobacillus ruminis Lactococcus lactis Odoribacter splanchnicus Parabacteroides distasonis Roseburia intestinalis Streptococcus spp. Veillonella dispar Veillonella parvula
Histidine degradation	Acidaminococcus intestini Alistipes finegoldii Alistipes shahii Bacteroides spp. Bifidobacterium longum Blautia coccoides Blautia producta Citrobacter youngae Clostridium sp. CAG-299 Enterobacter cloacae complex Escherichia coli Klebsiella spp. Streptococcus parasanguinis
Phenylacetate degradation	Citrobacter freundii Citrobacter youngae Enterobacter cloacae complex Escherichia coli Klebsiella spp.
Trans-cinnamate degradation	Escherichia coli Klebsiella spp.

In terms of the affected microbiome pathways containing vitamin B₁₂-dependent enzymes, deoxyribonucleotide biosynthesis, and de novo pyrimidine biosynthesis were both significantly negatively associated with neurodevelopmental metrics of the developmental skills domain (Table S7). Deoxyribonucleotide biosynthesis was the most significantly negatively associated KEGG module by t-score for cognition, while de novo pyrimidine biosynthesis was the most significantly negatively associated KEGG module by t-score for the IBQ-R effortful control and surgency scores and the most significantly positively associated KEGG module by t-score for the ICQ dull score (Table 2). Most of these effects were again found to be male-specific (Figure 3; Table S8). Several Enterobacteriaceae spp. associated with community type 2 were observed to contain the increased capacity for deoxyribonucleotide biosynthesis (Table 3). We further tested whether microbiome functions mediated the observed association between microbiome community type and cognition. The GSEA of the average casual mediation effects of the KOs revealed that both deoxyribonucleotide biosynthesis and pyrimidine deoxyribonucleotide biosyntheses partially mediated the association between microbiome community type and cognitive development (Table S9).

The KEGG modules methanogenesis, the methylaspartate cycle, and the ethylmalonyl pathway also were found to be significantly associated with several metrics in the developmental skills domain (Table S7). Methanogenesis was the most significantly positively associated KEGG module by t-score for surgency; the methylasparate cycle was the most significantly positively associated KEGG module by t-score for language skills and effortful control; and the ethylmalonyl cycle was the most significantly positively associated KEGG module by t-score with cortisol reactivity (Table 2). All effects on motor skills and difficult temperament for these microbiome function-neurodevelopmental associations were sex-dependent (Figure 3; Table S8). Blautia spp., the most associated microorganisms with community type 1, were found to be critical contributors for methanogenesis and methylaspartate cycle microbiome functions (Table 3). Eubacterium callanderi and Eubacterium limosum were also capable of methanogenesis, and interestingly, Enterobacteriaceae spp., the most associated microorganisms with community type 2, were additional contributors to the methylaspartate cycle and ethylmalonyl pathway (Table 3).

Gut microbiome antigens and pathogenicity signatures are broadly associated with both neurodevelopmental domains

Pro-inflammatory gram-negative microbial surface antigens were broadly negatively associated with both neurodevelopmental domains in male children. Specific anti-inflammatory, anti-microbial and competitive microbiome mechanisms were positively associated with both neurodevelopmental domains specifically in males.

Pathogenic Escherichia coli infection was the most significantly negatively associated KEGG pathway by t-score for the BSID-III motor composite score (Table 2), specifically in relation to fine motor skills (Table S7) and was also significantly associated with decreased BSID-III cognitive composite scores, and IBQ-R surgency and effortful control scores (Table S7). The eaeA gene (intimin virulence factor for EPEC/EHEC) was found to be present in the microbiome of > 20% of the male children but in none of the female children in this cohort. Therefore, associations between the neurodevelopmental domains and typical signatures of gram-negative, pathogenic bacteria were mostly specific to males (Figure 3; Table 2, Tables S7 and S8). Further, pathogenic E. coli infection was the most significant mediating pathway by t-score for the association between microbiome community type and BSID-III language composite score, in addition to significantly mediating the association between microbiome community type and cognition, both associations of which were specific to male children (Table S9). Moreover, when conducting GSEA of the t-scores from linear regression of the KO abundances between children with and without E. coli pathogenicity signatures (i.e., eaeA gene), vitamin B₁₂ biosynthesis was found to be the most significantly depleted KEGG module by t-score, suggesting that these opportunistic pathogens may exploit a vitamin B₁₂-deificent gut environment (Table S10).

Several microbiome-mediated mechanisms to control inflammation and microbial populations were positively associated with neurodevelopment in males (Figure 3; Tables S7 and S8). Secondary bile acid biosynthesis was found to be the most significantly positively associated KEGG pathway by t-score with cognition (Table 2), the most significant mediating KEGG pathway by t-score for the association between microbiome community type and cognition (Table S9), and the most significantly depleted KEGG pathway between children with and without E. coli pathogenicity signatures by t-score (Table S10). These functions were observed in relatively few bacterial spp.: Clostridium scindens, Dorea sp. D27, Eggerthella lenta, and Gordonibacter pamelaeae (Table 3).

Biosynthesis of dTDP-L-rhamanose is part of the O-Antigen nucleotide sugar biosynthesis KEGG pathway that was most significantly associated with fine motor skills by t-score (Table 2). Blautia spp., Bacteroides spp., Bifidobacterium longum, and several other Lachnospiraceae spp. were found to be capable of biosynthesizing dTDP-L-rhamanose (Table 3).

Finally, the KEGG pathway valine, leucine, and isoleucine degradation was significantly positively associated with language skills, specifically receptive communication, surgency, and effortful control, and was significantly negatively associated with the ICQ fussy/difficult, unadaptable, and unpredictable scores (Table S7). Blautia spp., Bacteroides spp., Bifidobacterium animalis and several other Lachnospiraceae spp. were all found to be capable of these functions (Table 3).

An association that affected both sexes was the abundance of the gram-positive microbial surface antigens, lipoarabinomannan (LAM) and teichoic acid, with developmental skills and temperament (Figure 3). LAM biosynthesis was the most significantly positively associated KEGG pathway by t-score with language skills, including both expressive and receptive communication, effortful control and surgency, and the most significantly negatively associated KEGG pathway by t-score with the IBQ-R negative affect and fussy/difficult temperament (Table 2). LAM biosynthesis was also significantly positively associated with cognitive and motor development, and significantly negatively associated with unadaptable and unpredictable temperaments, making it the most broadly significantly associated microbial metabolic pathway across both neurodevelopmental domains (Table S7). Although most of these associations were sex-independent, the effect on language skills was significantly stronger in males, and the impact on negative affect and unpredictable temperament was significantly stronger in females (Figure 3; Table S8). LAM biosynthesis capabilities were present in Corynebacterium spp. and Bifidobacterium animalis (Table 3). Additionally, teichoic acid biosynthesis was the most significantly positively associated KEGG pathway by t-score with gross motor skills (Table 2). Only the temperament impacts of teichoic acid were exclusive to males, whereas its associations with language and motor skills were sex independent (Figure 3; Tables S7 and S8). Blautia coccoides, Lactobacillus spp., Lactococcus spp. and Streptococcus spp. all possessed genes for teichoic acid biosynthesis (Table 3).

Purine, histidine, and phenylalanine metabolism modulated by the gut microbiome are associated with difficult temperament

In contrast to the developmental skills domain, there was no singular microbiome function associated with all aspects of the difficult temperament domain. The closest microbiome function to being key to difficult temperament was purine degradation, which was significantly associated with 4/5 of its facets (Table S7). Purine degradation was the most significantly negatively associated KEGG module by t-score with the IBQ-R negative affect score and ICQ unadaptable score (Table 2). In contrast to the previous findings in the developmental skills domain, most of the significant associations between purine degradation and difficult temperament metrics were exclusive to females (Figure 3; Table S8). This microbiome function was again found to be contributed by Blautia spp. (Table 3).

Other metabolic pathways related to difficult temperament were histidine metabolism and phenylalanine metabolism (Table S7). The significant associations observed between neurodevelopmental outcomes and histidine metabolism were mostly sex-independent, whereas phenylalanine metabolism was mainly male-specific in a complex fashion (Figure 3; Table S8). Histidine metabolism was the most significantly positively associated KEGG pathway by t-score for the ICQ fussy/difficult and unpredictable scores, and the histidine biosynthesis KEGG module specifically was the most significantly positively associated KEGG pathway by t-score for unpredictable temperament (Table 2). Histidine metabolism is contributed by multiple microorganisms, although Enterobacteriaceae spp. are relevant to our findings (Table 3). Phenylacetate degradation was the most significantly negatively associated KEGG module by t-score for the BSID-III motor composite and expressive communication scaled scores (Table 2), but also was significantly negatively associated with negative affect and unpredictable temperament while being positively associated with the ICQ dull score (Tables S7 and S8). Trans-cinnamate degradation was the most significantly positively associated KEGG module by t-score for cortisol reactivity (Table 2), but also was significantly negatively associated with fussy/difficult and unpredictable temperaments (Table S7). Enterobacteriaceae spp. were the exclusive contributor of these phenylalanine-derived metabolic pathways (Table 3).

Metabolomics analysis of fecal samples corroborates a role for secondary bile acids in the cognitive outcomes of male children

Targeted mass spectrometry-based metabolomics analysis of the infant fecal samples was conducted to validate the observed associations between key microbiome functional KEGG pathways and neurodevelopmental outcomes via shotgun metagenomics. The measured metabolites included vitamin B₁₂ (cyanocobalamin & hydroxocobalamin), methylaspartate cycle metabolites (acetate & glutamate), ethylmalonyl pathway metabolites (propionate & succinate), 49 bile acids (primary, secondary, and glyco/tauro-conjugated subclasses), branched-chain amino acids (valine, leucine & isoleucine), branched-chain fatty acids (isobutyrate, isovalerate & 2-methylbutyrate), purine degradation metabolites (guanine, adenosine, inosine, hypoxanthine, xanthine, uric acid & urea), histidine metabolites (histidine, histamine, urocanic acid, imidazoleacetic acid & imidazoleproprionic acid), and phenylalanine, phenylacetate & trans-cinnamate.

Of these metabolites, vitamin B₁₂, most purine degradation metabolites (except for inosine & urea), most histidine metabolites (except for histidine & urocanic acid), phenylalanine & trans-cinnamate were below the limit of detection across all samples (Table S11). From regression analysis that adjusted for gestational age at birth, birthweight, and age at assessment visit, while also including sex as a moderator, acetate was found to significantly negatively associate (p < 0.05), whereas glutamate was found to significantly positively associate in males, with the BSID-III motor composite score (Table S12). Glutamate also overall significantly negatively associated with IBQ-R negative affect (Table S12). Additionally, propionate was significantly negatively associated and succinate significantly positively associated with the BSID-III fine motor scaled score in males (Table S12). Succinate was positively associated with IBQ-R surgency in a male-specific manner (Table S12). For purported temperament-associated metabolites, urocanic acid was significantly positively associated with IBQ-R surgency and effortful control with male-specificity (Table S12). Finally, phenylacetate was significantly positively correlated with the BSID-III cognitive composite score (Table S12).

Of all metabolite classes measured in fecal samples, bile acids had the highest number of significant associations with infant neurodevelopmental outcomes. In particular, 17 distinct bile acids significantly negatively correlated with cognition in male children (Figure 4; Table S12). There were also bile acids that were significantly associated with language & motor skills, surgency & effortful control, and observed stress regulation (i.e., cortisol reactivity after a stressor), indicating that bile acids were associated with all metrics of the developmental skills domain.

Figure 4. The significant (p < 0.05) estimated marginal linear trends (emtrend) between fecal bile acid content and the Bayley Scales of Infant Development III cognitive composite score by sex in a cohort of Black children living in low socioeconomic and high stress environments (average age 7.6 months old; n = 28), after adjustment for gestational age at birth, birthweight, and age at assessment visit.

Cholic acid-7-sulfate, lithocholic acid-3-sulfate and tauroursodeoxycholic acid all significantly negatively associated with the BSID-III language composite score in males (Table S12). Glycolithocholic acid and glycodehydrocholic acid were significantly negatively associated with fine motor and gross motor skills in male children, respectively, whereas glycodeoxycholic acid was significantly negatively correlated with fine motor skills and glycoursodeoxycholic acid was significantly positively correlated with gross motor skills with sex-dependence (Table S12). Cholenic acid was negatively correlated with both surgency and effortful control in males, and glycocholic acid, taurocholic acid-3-sulfate, taurohyodeoxycholic acid and 7-oxodeoxycholic acid additionally all were negatively associated with surgency in male children (Table S12). Regarding cortisol reactivity, chenodeoxycholic acid-3-sulfate and glycodehydrocholic acid were significantly negatively correlated with cortisol reactivity in males (Table S12).

Finally, there were bile acids that significantly associated with aspects of difficult temperament. Hyodeoxycholic acid and taurocholic acid-3-sulfate were significantly negatively correlated with negative affect in male children, oxocholic acid was negatively associated and taurolithocholic acid was positively associated with negative affect without sex-dependence, and chenodeoxycholic acid was overall negatively correlated with ICQ fussy/difficult temperament (Table S12).

Discussion

This study was the first to examine concurrent relationships between gut microbiome composition and function and multiple dimensions of neurodevelopment in a sample of Black infants (n = 28; 1 sample per participant) living in high stress, under-resourced, urban environments. Our results indicate that the biosynthesis of vitamin B₁₂ is a key microbiome function associated with all measured facets of neurodevelopment including cognition, language skills, motor development, surgency, effortful control, and observed stress regulation (i.e., cortisol reactivity after a stressor). Microbial pathways containing vitamin B₁₂-dependent enzymes were additionally affected. These significant findings were found in male children only, and indeed > 20% of male children in this cohort contained E. coli with pathogenic signatures, necessitating the significantly associated anti-microbial, anti-inflammatory, and resource competition mechanisms. Blautia spp. were identified as a top contributor of microbiome functions found to be important for neurodevelopment.

Vitamin B₁₂ is an essential micronutrient in humans that plays a critical role in brain development both in and ex utero, including processes of cell division and growth, myelination, synaptogenesis and neurotransmitter synthesis.¹⁷ We cannot definitively rule-out vitamin B₁₂ deficiency in this cohort as blood samples were not collected; however, vitamin B₁₂ deficiency is rare in the United States, especially for people of African ancestry who have a genetic predisposition to higher absorption of vitamin B₁₂, thus yielding significantly higher serum levels than other ancestries.^18,19 The gut microbiome is not thought of as a major contributor to host vitamin B₁₂ levels, as vitamin B₁₂ is absorbed in the small and not large intestine.^20,21 However, host vitamin B₁₂ intake could potentially impact the gut microbiome, as dietary vitamin B₁₂ is inefficiently absorbed in the small intestine, normally leading to sufficient quantities in the colon.²¹

Measuring the total colonic contents of vitamin B₁₂ is challenging, as > 98% of ingested cobalamin is converted into cobalamin analogues or corrinoids.^20,22 Using a specially developed column purification method combined with liquid chromatography-mass spectrometry, Allen and Staber determined that the average concentration of cobalamin in adult feces was 1309 ng/g,²² which is below the limit of detection for standard mass spectrometry-based assays. Gut microorganisms can transform cobalamin into at least eight distinct corrinoids through a process termed corrinoid remodeling, which presents a competitive advantage, because each corrinoid can only be used by a fraction of the species in the gut microbiota since corrinoid-dependent microbial enzymes have native specificity for their preferred cofactors.²⁰ Therefore, although we did not detect vitamin B₁₂ in any of the infant fecal samples collected for this study, this does not negate the results from the shotgun metagenomics analysis. Rather the lack of detection is likely indicative of a limitation of metabolomics analysis, which measures only the net metabolomic output and not the total input (i.e., ingested vitamin B₁₂) or total output (i.e., vitamin B₁₂ produced by microorganisms).

Vitamin B₁₂ is a necessary cofactor for enzymes essential to the viability of most intestinal bacteria, including ribonucleotide reductase (nrdJ – essential for deoxyribonucleotide biosynthesis).^20,23 However, many intestinal bacteria have also evolved to contain vitamin B₁₂-independent enzymatic equivalents, including nrdAB, nrdEF and nrdD, and thus auxotrophy for vitamin B₁₂ is a rare but not unfeasible trait of gut-residing bacteria.²³ Even in the absence of auxotrophy, enzymatic flexibility in these essential functions may nonetheless be a selection factor for bacteria in a vitamin B₁₂-deficient gut. Enterobacteriaceae spp. encode a large diversity of ribonucleotide reductases which allow them to adapt to various environments, and a knock-out of just one of its ribonucleotide reductases reduced the ability of adherent-invasive Escherichia coli to colonize the gastrointestinal tract of mouse models and express virulence genes, demonstrating not only a loss of competitive advantage but a change in microbial behavior through altered expression of ribonucleotide reductases.^24,25 As part of these changes in behavior for gut microorganisms, vitamin B₁₂ is also a necessary cofactor for multiple nonessential enzymes catalyzing methyl transfer reactions. These enzymes include trimethylamine methyltransferase for methanogenesis from trimethylamine, methylaspartate mutase for consumption of glutamate and acetyl-CoA as part of the methylaspartate cycle for anabolic processes, and methylmalonyl-CoA mutase for bi-directional interconversion of propanoyl-CoA and succinyl-CoA to control CO₂ and propionate production via the ethylmalonyl pathway.²⁰ The net result would potentially indicate altered trimethylamine, CH₄, glutamate, acetate, propionate, succinate, and CO₂ levels in the gut, with trimethylamine, glutamate and propionate being microbial metabolites particularly associated with neurodevelopmental/neurological disorders in several studies.^26–29

Due to their volatile nature, we could not measure CO₂, CH₄ and trimethylamine by standard mass spectrometry-based metabolomics analysis in fecal samples. For fecal concentrations of glutamate, acetate, propionate, and succinate, however, significant associations were only detected with motor skills, surgency, and negative affect. The relative lack of significant associations for these metabolites and neurodevelopment may be due to some of the inherent limitations of mass spectrometry-based metabolomics analysis, including host interference via metabolite consumption by intestinal epithelial cells and absorption into the systemic circulation, daily metabolite fluctuations from dietary changes, and the aforementioned restriction of measurement to net metabolic output.^30,31 Alternatively, vitamin B₁₂-dependent microbiome metabolic pathways may impact the developing microbiome more from potentially fostering pathogenicity than through altering metabolic output. Without blood samples, we were unable to evaluate systemic immunity, the other main mechanism by which the microbiome interacts with its host. It is probable that the Enterobacteriaceae-dominant community type putatively enabled by vitamin B₁₂-deficiency promoted a chronic low-level systemic pro-inflammatory host response, which has been found to hinder neurodevelopment in both animal models^5,32 and human observational studies^33–35. This discovered mechanistic link between vitamin B₁₂-deficiency and pathobiont colonization leading to inflammation and suboptimal neurodevelopment should be examined in future work.

Current dietary guidelines have mainly focused on preventing nutrition deficiencies in the host, but there is increasing evidence that the microbiome itself has nutritional needs. Variations in dietary protein, fat, phytochemicals, food additives/artificial sweeteners and micronutrients have all been associated with compositional and/or functional deviations of the microbiome,^36,37 and understanding the appropriate dietary quantities for microbiome nutrition may be especially critical during early life to both optimize microbiome successional, and thus host developmental, trajectories. As an example in the context of our study, among racial and ethnic groups in the United States, Black women have the lowest vitamin B₁₂ dietary intake.³⁸ Although efficient intestinal absorption may preclude host vitamin B₁₂ deficiency, it may also result in a relatively low amount of vitamin B₁₂ reaching the colonic microbiome potentially impacting the maternal microbiome that is vertically transmitted from mother to infant. Indeed, Selma-Royo and colleagues demonstrated that maternal vitamin B₁₂ intake was positively correlated with Bifidobacterium and several Lachnospiraceae spp. in infant stool shortly after birth.³⁹ Thus, the interplay between host nutrition and microbiome succession is a topic deserving of future research, and vitamin B₁₂ may particularly serve as a scalable nutrition intervention to protect early neurodevelopment, especially for Black male infants living in under-resourced communities.

We observed that deviations in microbiome structure were associated with neurodevelopmental outcomes in male and not in female infants. This finding is in alignment with another study by Tamana et al. that examined the relationship between fecal microbiome composition and BSID-III scores in a predominantly White sample of 1–2-year-old children from the Canadian Healthy Infant Longitudinal Development cohort, which determined that a Bacteroidota-dominant community type was associated with better cognitive and language composite scores exclusively in male children.⁴⁰ Males are known to be at a significantly higher risk for developmental disorders, learning disabilities and cognitive impairment than females.^41,42 It is of note that that results from both the present study and the study conducted by Tamana et al. did not yield differences in the microbiome between the sexes, but instead demonstrated sex-specific responses to microbiome dysbiosis (Figure 2a). Results from two other observational studies^43,44 showed that early microbiome composition (<3 months old) was prospectively associated with maternal report of temperament at 6–12 months of age, underscoring the significance of the microbiome in early life to brain development. Notably, in all of these studies hypothesis testing was limited to microbiome composition. Our study thus extends the extent research by identifying specific microbiome functions associated more with neurodevelopment.

Differences in the early caregiving environment of male and female infants may partially explain the differential association between microbiome and neurodevelopment. For example, preliminary data from our study indicate an earlier introduction of solids for males than females that occurred at about 4 months of age. This finding is consistent with results from other studies showing that male children are significantly more likely to have an early introduction of solid foods than female children, usually because they are perceived as “hungrier” and “less satisfied” with breastfeeding alone.^45,46 Solid food introduction increases both the abundance and metabolic output of the gut microbiome due to the increase in fermentable carbohydrates.⁴⁷ Therefore, an imbalance of gram-negative bacteria might become particularly inflammatory once their numbers and thus total antigenic content in the gut increases, and an imbalance of metabolites might begin to impact the host at increased concentrations. Further, nutritional needs of the microbiome, including vitamin B₁₂, may change to compensate for the elevated bacterial load. Differding and colleagues have demonstrated that the early introduction of solid foods (≤3 months old) altered the abundance of gut microbial taxa and fecal SCFA concentrations at 3 and 12 months old when compared to controls,⁴⁸ and studies have associated the early introduction of solid foods with obesity, immune-mediated conditions (e.g., allergy) and increased risk of respiratory and GI infections.⁴⁹ Indeed, > 20% of male children in our cohort contained pathogenicity markers for EPEC/EHEC compared to none of the female children. Therefore, it is of interest in future work to examine whether sex-dependent feeding attitudes and practices intersect with the microbiome to influence developmental outcomes.

Fortunately, several microbiome mechanisms that bolster innate immune defenses against opportunistic pathogens were uncovered in this study that have potential to improve neurodevelopmental outcomes. Hydrophobic bile acids produced from 7α-dehydroxylation or 3α-dehydrogenation have been found to be potent antimicrobials particularly against gram-negative bacteria.^50,51 Metabolomics analysis of infant fecal samples confirmed that at least 17 different bile acids were correlated with cognitive outcomes in males in this study (Figure 4). Of these 17 bile acids, seven conjugated bile acids and cholic acid (a primary de-conjugated bile acid) were negatively correlated, which could indicate a buildup of secondary bile acid biosynthesis reactants in males with suboptimal cognitive outcomes due to a lack of microbial activity. Additionally, the bacterial surface sugar dTDP-L-rhamanose has been found to attenuate LPS-mediated inflammation and induce the production of anti-inflammatory cytokines.^52–54 Lastly, valine, leucine, and isoleucine degradation has been found to deplete resource pools which modulate virulence gene expression of several bacterial pathobionts,⁵⁵ although fecal metabolomics analysis could not corroborate this association.

However, these findings do not preclude specific microbiome mechanisms from benefiting female children, particularly in terms of temperament outcomes, which include purine degradation, histidine metabolism, and phenylalanine metabolism. Purines are degraded to uric acid, of which serum concentrations that are both too low and too high have been associated with cognitive impairments, but the beneficial effects of uric acid are attributed to its antioxidant properties.⁵⁶ Histidine is a psychoactive amino acid acting as a precursor to the histamine neurotransmitter^57,58 and imbalances in phenylalanine-derived metabolites have been previously correlated with ASD.^59–62 The studied purine-, histidine- and phenylalanine-derived metabolites were largely below the limit of detection in infant fecal samples by standard mass spectrometry-based metabolomics analysis, which may indicate that these compounds were subject to previously described host interference or further breakdown by microbial metabolism. A final important microbiome mechanism to note are the gram-positive surface structures LAM and teichoic acid that were broadly associated with developmental skills and temperament for both male and female children. These antigens are toll-like receptor (TLR) 2-agonists that have been demonstrated to promote colonization through immunosuppressive effects.^63–65

We identified Blautia spp. as the most dominant microbial genus for community type 1 that was positively associated with cognition and language skills, and also contributed many of the identified key microbiome functions: vitamin B₁₂ biosynthesis, regulation of glutamate and acetate via the methylaspartate cycle, methanogenesis from trimethylamine, the anti-inflammatory microbial antigens teichoic acid and dTDP-rhamanose, degradation of branched-chain fatty acids to deplete resource pools for incoming potential pathogens, and purine degradation to the antioxidant uric acid. Blautia has been previously found to be depleted in patients with ASD,⁶⁶ and Sen and colleagues demonstrated that administration of Blautia stercoris attenuated social, repetitive, and anxiety-like behavior in a mouse model.⁶⁷ Reduced Blautia abundance has also been associated with major depressive disorder⁶⁸ and cognitive impairment⁶⁹ in adults. Interestingly, Blautia has been implicated in vertical transmission from mother to child, and breastfeeding and maternal weight may impact the success of transmission.^70,71 Taken together, Blautia spp. are a possible therapeutic target for improving neurodevelopmental outcomes in children who may be at risk for suboptimal development due to familial or environmental factors.

The observed association between gut microbiome dysbiosis and neurodevelopmental outcomes in this cohort of Black American children living in high-stress, resource constrained environments likely was caused by multiple factors. These factors include mode of delivery, formula amount/type, amount/type of the solid food, antibiotic usage and infections, other medications and infant health determinants (e.g., growth), exposure to environmental toxins (e.g., secondhand smoke, pollutants), psychosocial stress, and health of the mother, other caregivers and siblings that could affect the passage of microbes to the infant.^8–10 The goal of the present study was to examine if and how dysbiosis correlated with neurodevelopment, in order to determine which microbes and their functions could have mechanistic impact and thus therapeutic utility. It would be of interest in future work, however, to examine the relative influence of these factors on the development of dysbiosis for children, which could yield further insights into preventative strategies to mitigate adverse microbiome successional trajectories.

We note that scores for the infants in our cohort were in the average range for neurodevelopment, but with sufficient variability to observe associations with the microbiome. Testing associations within a cohort of Black infants living in environments that are comparable with regard to stress, financial strain, and other environmental injustices, allows the study of both risk and resilience. Understanding how the infant microbiome is natively optimized to promote development in the context of environment stressors is equally important to developing scalable preventive interventions that are relevant for youth living under such circumstances. Further investigation of vitamin B₁₂ as a nutritional target and Blautia spp. as a candidate live biotherapeutic are examples of potential interventions to support neurodevelopment in this population of children. Future work should address this study’s limitations by including larger cohorts with greater diversity in race and ethnicity and multiple geographic areas that vary in urbanicity to extend the generalizability of our findings. Additionally, inclusion of larger cohorts with greater diversity in stress exposure would provide a more nuanced approach to the measurement of race-related stressors. For example, including Black families living in higher resourced communities would allow testing for the impact of subsistence-related stress (e.g., housing stress) versus interpersonal discrimination stress on the associations between the gut microbiome and neurodevelopment. Further, longer-term follow-up should be conducted to determine if and how associations observed at earlier ages impact later school readiness and early childhood functioning and well-being.

Participants and methods

Participants

Participants were drawn from the Nutrition and Pregnancy Study (NAPS: NCT02647723), an ongoing double-blind randomized controlled trial of fatty acid supplementation during pregnancy in Black women with Medicaid insurance. The NAPS study is aimed at reducing health disparities by improving maternal stress regulation during pregnancy, by testing the effect of omega-3 fatty acid supplementation among women reporting low levels of sea (i.e., saltwater or ocean-dwelling) fish consumption. Inclusion criteria included self-identification as Black, age between 18–34 years, receipt of public assistance (i.e., Medicaid insurance), and less than two servings of sea fish per week. Medicaid is a joint state and federal government health insurance program in the United States for individuals and families with low-incomes.¹⁴ Exclusion criteria comprised serious medical complications, regular use of steroid medications, blood thinners or anti-coagulants, psychotropic medications, substance use, and allergy to fish, soy, strawberries or iodine. Participants eligible for enrollment were identified through medical records at the University of Chicago and the Friend Family Health Center, a federally qualified health center affiliated with the University of Chicago. Approval for all study procedures was obtained from the University of Chicago Institutional Review Board (IRB# 150392) which follows the Declaration of Helsinki. Written informed consent was obtained prior to data and sample collection. All infants from the NAPS study who were age-eligible (<1.5 years old) at the time of ethics approval for fecal specimen collection were approached for consent from the parent. Of the parents approached, 71% consented to child fecal specimen collection, and none of the data from these children were excluded for this study.

The Great Cities Community Area Hardship Index^15,16 was calculated for participant families based upon their zip code at the time of study recruitment. Data from the 2016 U.S. Census Bureau’s American Community Survey were used to calculate index values by census tracts, and Geographic Information Systems (GIS) software was used to aggregate census tract data and geographies into Chicago Community Area boundaries. The economic hardship score is an average of six variables: 1) unemployment (over the age of 16 years); 2) education (over 25 years of age without a high school diploma); 3) per capita income level; 4) poverty (households below the poverty level); 5) crowded housing (housing units with more than one person per room); and 6) dependency (population under 18 or over 64 years of age), which are standardized on a scale from 0 to 100.

Covariate data and sample collection

Covariate data was collected from medical records or by maternal report during the study visits when the infant was four months old and nine months old, which included infant gestational age at birth, birthweight, date of birth, sex, delivery mode, breastmilk or formula feeding, receipt of solid foods (i.e., cereal) in the past week, amount of sleep at night (19:00 to 07:00) in hours, amount of sleep during the day (07:00 to 19:00) in hours, number of night wakenings, and length of time to put baby to sleep at night in minutes. Fecal samples were collected from the children at each neurodevelopmental assessment visit (4 months, 9 months, and 12-months old), and a single fecal sample from each study participant with the corresponding neurodevelopmental evaluation results was utilized for analysis. The latest visit available for each participant was chosen for analysis, and a total of 28 children provided fecal samples for this study.

Developmental skills and temperament assessment

Neurodevelopment was assessed in infancy (<1.5 years old) and included the following: Bayley Scales of Infant Development III (BSID-III) for cognition, language (receptive and expressive communication), and motor (fine and gross) skills; maternal report on the Infant Behavior Questionnaire Revised (IBQ-R) for surgency (i.e., extroversion), negative affect and effortful control (i.e., orienting attention/regulation); maternal report on the Infant Characteristics Questionnaire (ICQ) for fussy/difficult, adaptability, predictability and dullness; and infant cortisol reactivity to the Face-to-Face Still-Face paradigm. For the BSID-III, the composite cognition, language and motor scores were computed for analysis, and the receptive communication, expressive communication, fine and gross motor scaled scores also were used.⁷² The BSID-III has shown good concurrent validity among Black infants; predictive validity is less optimal for infants of all races.^73–75 The ICQ contains 32 seven-point items that assess temperamental characteristics, including fussy-difficult, unadaptable, unpredictable, and dull characteristics.⁷⁶ In samples of families experiencing economic strain, parental report of temperament for Black and White infants using the ICQ showed strong internal consistence and good predictive validity.⁷⁷ The IBQ-R very short form contains 37 items that load on three dimensions: Surgency/Extraversion, Negative Affectivity, and Orienting/Regulation (i.e., Effortful Control);⁷⁸ good internal consistency is reported for infants of different races and from different economic environments.⁷⁸ For all metrics in the BSID-III, IBQ-R and ICQ, the stated scores were used directly as numeric dependent variables in multiple regression analysis, instead of stratifying the study participants into optimal and suboptimal developmental groupings.

Infant observed stress regulation

Infant response to the Face-to-Face Still-Face paradigm (FFSF) was used to measure infant observed stress regulation.⁷⁹ The FFSF is a standard laboratory procedure comprising three 2-minute episodes: (1) mother playing typically with her seated infant; (2) mother maintaining a neutral expression with no vocalization; and (3) mother returning to typical play. Mothers were asked that their infants should not consume milk products (human or animal) for 1 h prior to and during saliva collection. Saliva was collected pre-FFSF, and 20, 40, and 50-min post-FFSF by swabbing each infant’s mouth with an unflavored dental roll for several minutes. Samples with sufficient saliva, after centrifuging the dental rolls at 3,000 rpm for 10 min, were assayed in duplicate using the Salimetrics HS Salivary Cortisol EIA Kit for unbound cortisol (Salimetrics, 1–3002). Cortisol reactivity was determined as the area under the curve with respect to increase (AUCI) calculated by the composite trapezoid rule via R package DescTools version 0.99.44 of the four cortisol measurements after subtraction of the pre-stressor value.

Illumina 16S rRNA gene sequencing and processing

Participant fecal samples were submitted to the Environmental Sample Preparation and Sequencing Facility at Argonne National Laboratory for genomic DNA extraction and Illumina 16S rRNA gene sequencing.^80,81 Data retrieved from the facility (NCBI SRA: PRJNA895372)⁸² were subsequently processed and merged by the sample inference tool DADA2⁸³ from within QIIME2 version 2019.10.⁸⁴ The naïve Bayes classifier built on the SILVA database⁸⁵ release 138 implemented in QIIME2 was utilized for classification to the genus level. The α-diversity metrics of richness and Shannon diversity were computed by R package iNEXT⁸⁶ version 2.0.20. For β-diversity analysis, the taxonomic levels of phylum, class, family, and genus were individually considered. Community typing was conducted by the probabilistic method of Dirichlet Multinomial Mixtures⁸⁷ for each of the taxonomic levels with 1 to 27 (n-1) community types considered and the optimal solution selected by the minimum Laplace value using R package DirichletMultinomial version 1.36.0.

Shotgun metagenomics sequencing and processing

Participant fecal samples were submitted to the Duchossois Family Institute at the University of Chicago for genomic DNA extraction (QIAamp PowerFecal Pro DNA Kit, Qiagen 51804), library generation (QIAseq FX Library Kit, Qiagen 180479), and shotgun metagenomics sequencing on the Illumina NovaSeq 6000 platform using the 2 × 150 Paired End read cassette. Data retrieved from the facility (NCBI SRA: PRJNA895372)⁸² were subsequently processed and merged by the bioBakery whole metagenome shotgun workflow version 3.0.0.⁸⁸ Briefly, this workflow utilizes the KneadData tool to perform quality control and remove host reads, the MetaPhIAn tool to generate taxonomic abundance profiles, and the HUMAnN tool to obtain microbial gene abundance profiles stratified by contributing microbial taxa. The MetaPhIAn tool uses the National Center for Biotechnology Information (NCBI) genome catalog⁸⁹ and taxonomic database,⁹⁰ and the HUMAnN tool uses the Kyoto Encyclopedia of Genes and Genomes (KEGG) database.⁹¹ The relative abundance data produced by this pipeline was center-log ratio transformed with zeroes imputed by the nonparametric multiplicative simple method through R package zCompositions version 1.4.0.1, to allow standard statistical testing to be conducted on the inherently compositional data.

Metabolomics analysis

Metabolomics analysis was conducted at the Duchossois Family Institute (DFI) at the University of Chicago. Compounds were specifically chosen for study in order to validate the microbiome functional KEGG pathways found to significantly correlate with neurodevelopmental outcomes from shotgun metagenomics analysis of infant fecal samples. The infant fecal metabolome was analyzed across three mass spectrometry platforms to capture qualitative levels of gut-derived metabolites with varying physiochemical properties such as hydrophobicity, size, and charge. In brief, metabolites were extracted with organic solvent, dried down and resuspended for direct analyses or derivatization. Previously, all compounds have been validated by the DFI through retention time and fragmentation comparison to standards and available databases. Gas chromatography-mass spectrometry (GC-MS) was used to detect compounds following derivatization with pentafluorobenzyl bromide (PFBBr)⁹² and trimethylsilyl-methoxamine (TMS-MOX)^93,94 in two separate reactions. Acetate, glutamate, propionate, succinate, valine, leucine, isoleucine, isobutyrate, isovalerate and 2-methylbutyrate were qualitatively analyzed by normalized peak area following PFBBr derivatization and detection by negative collision induced-GC-MS ((-)CI-GC-MS, Agilent, 8890). Positive ion electron impact-GC-MS ((+)EI-GC-MS, Agilent, 7890B) was used to detect guanine, adenosine, inosine, hypoxanthine, xanthine, uric acid, urea, histidine, histamine, urocanic acid, imidazoleacetic acid, imidazole propionic acid, phenylalanine, phenylacetate, trans-cinnamate, cyanocobalamin and hydroxocobalamin following TMS-MOX derivatization. With the use of negative mode liquid chromatography-electrospray ionization-quadrupole time-of-flight-MS ((-)LC-ESI-QTOF-MS, Agilent, 6546), 49 bile acids from the primary, secondary and glyco/tauro-conjugated subclasses were analyzed.⁹⁵ In addition to retention time validation, the standard intact and fragment masses are routinely detected with differences < 5 ppm compared to calculated values. Metabolite normalized peak abundances are recorded in Table S11.

Statistical analysis

All analysis was conducted in R statistical software version 4.1.2, and plots were generated using R packages ggplot2 version 3.3.5, corrplot version 0.92, interactions version 1.1.5 and patchwork version 1.1.3. For each neurodevelopmental outcome measurement, multiple linear regression models were constructed to determine if richness, Shannon diversity, community type, KEGG orthologies (KOs) and metabolites were significantly related. Adjusted cofactors for all models included gestational age at birth (weeks), birthweight (grams), age at evaluation visit (days) and sex. For the cortisol reactivity outcome, time of day (hours passed midnight) was included as an additional covariate. To determine how sex moderated the association between the microbiome variable and the neurodevelopmental outcome, the interaction between the microbiome variable and sex was included as an additional model term. All model variables were standardized using R package standardize version 0.2.2 prior to fitting the models. Model assumptions were assessed by the Shapiro-Wilk test of normality on the model residuals, the Rainbow test for linearity and Harrison-McCabe test for heteroscedasticity from R package lmtest version 0.9.39, and the Bonferroni outlier test and variance inflation factors to assess collinearity (maximum <5 considered not collinear) from R package car version 3.0.12. Models that both met the statistical assumptions and were overall significant (p < .05) by ANOVA (Pr > F) were considered. If richness, Shannon diversity, community type and metabolites were found to be significantly (p < .05) related to a neurodevelopmental outcome measure, the term significance (Pr > |t|), coefficient and associated 95% confidence interval, and model adjusted R² are reported. Results and assumption tests of all models (significant or non-significant) are included in Table S6 (α-diversity and β-diversity) and Table S12 (metabolites). For the KOs, the organized t-scores were utilized for gene-set enrichment analysis (GSEA) via R package clusterProfiler version 4.4.4 with a minimum gene-set size of 3 for KEGG pathways and 2 for KEGG modules. If a KEGG module or pathway was found to be significantly (p < .05; false-positive rate < 1%) related to a neurodevelopmental outcome measure, the normalized enrichment score and p value are reported. GSEA results are included in Table S7 (overall) and Table S8 (sex-moderated).

When multiple linear regression model assumptions failed to be met, the equivalent log-linked Gamma regression models were considered. For these models, assumptions were evaluated using R package DHARMa version 0.4.5 plus the calculation of variance inflation factors to assess collinearity as done previously. Overall model significance was assessed by the chi-square test against the null model (Pr > Chi), and models that both met the statistical assumptions and were overall significant (p < .05) were considered. If richness, Shannon diversity, community type and metabolites were found to be significantly (p < .05) related to a neurodevelopmental outcome measure, the term significance (Pr > |t|), coefficient and associated 95% confidence interval, and Nagelkerke’s R² calculated by R package fmsb version 0.7.3 are reported. Results and assumption tests of all models (significant or non-significant) are also included in Table S6 (α-diversity and β-diversity) and Table S12 (metabolites). GSEA from the KO data was additionally conducted as above and results are included in Table S7 (overall) and Table S8 (sex-moderated).

Sex differences in the microbiome were evaluated as above, with linear or log-linked Gamma regression utilized for richness and Shannon diversity, and binomial regression utilized for microbiome community type. The methodology for binomial regression was the same as for log-linked Gamma regression, and results of all models (significant or non-significant) are included in Table S6. GSEA from the KO data was additionally conducted as above using simple linear regression for the microbiome community types (Table S9) and between children with and without E. coli pathogenicity signatures (Table S10), defined as the presence or absence of the eaeA intimin virulence factor for EPEC/EHEC (K12790). Differences in the neurodevelopmental outcome measures by sex and other covariates were also determined using simple linear or log-linked Gamma regression models in a similar fashion as above, with results of all models (significant or non-significant) included in Table S2. Sex differences for clinical covariates were ascertained using simple linear, log-linked Gamma or binomial regression as above, with results of all models (significant or non-significant) included in Table S3. Finally, differences in maternal DHA levels at 36 weeks gestation between the randomized controlled trial (NAPS: NCT02647723) groupings of prenatal DHA supplementation and placebo were determined through simple linear and log-linked Gamma regression as above, with results of all models (significant or non-significant) included in Table S1.

Mediation analysis for the associations between microbiome community type and BSID-III cognitive or language composite score for each KO was determined using the quasi-Bayesian approximation via R package mediation version 4.5.0. Linear models that adjusted for gestational age at birth, birthweight, age at assessment visit and sex were utilized, except for when the BSID-III language composite score was the outcome due to its prior non-linear association with microbiome community type, in which case the equivalent log-linked Gamma model was used. The resulting average casual mediation effects (ACME) were ranked and inputted into GSEA as done previously.

The partial Pearson correlation of neurodevelopmental outcome measures after adjustment for gestational age at birth, birthweight, age at evaluation visit and sex was conducted using R package ppcor version 1.1 and is presented as a heatmap ordered by hierarchical clustering via the ward linkage method (Figure 1). Permutation around medoids (i.e., k-medoids) clustering of the partial Pearson correlation matrix was conducted using R package fpc version 2.2.9 with 1 to 16 (n-1) clusters considered and the optimal solution selected by the maximum average silhouette width.

Acknowledgments

We acknowledge the Nutrition and Pregnancy families for their partnership, collaboration, and generosity. We thank the Duchossois Family Institute (DFI) for providing their expert technical services. We thank the Center for the Science of Early Trajectories (SET) for providing research coordinator resources.

Disclosure statement

Bree Andrews, MD/MPH is an equity partner in Preeme+You, a social benefit corporation, that support using mobile technology in the NICU to improve parent engagement and physician communication. There is no discussion of mobile technology in this manuscript. Erika Claud, MD has served as an expert witness for legal proceedings associated with the outcomes of infants in the neonatal intensive care unit unrelated to this research. The remaining authors have no competing interests they wish to declare.

Data availability statement

The data that support the findings of this study are openly available in NCBI SRA at http://www.ncbi.nlm.nih.gov/bioproject/895372, reference number PRJNA895372, in addition to the provided supplementary tables.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/19490976.2023.2298697

Additional information

Funding

This work was supported by the US National Institutes of Health under Grants R01 HD084586 and R01 HD105234; and the University of Chicago Comer Development Board Grant.