Neurodevelopment correlates with gut microbiota in a cross-sectional analysis of children at 3 years of age in rural China

We investigated cross-sectional associations between children’s neurodevelopment and their gut microbiota composition. Study children (36 months of age) lived in rural China (n = 46). Neurodevelopment was assessed using the Bayley Scales of Infant Development, 2nd Edition, yielding the Mental Developmental Index (MDI) and Psychomotor Developmental Index (PDI). Children's gut microbiota was assessed using 16S rRNA gene profiling. Microbial diversity was characterized using alpha diversity patterns. Additionally, 3 coabundance factors were determined for the 25 most abundant taxa. Multivariable linear regression models were constructed to examine the relationships between Bayley scores (MDI and PDI) and children's gut microbiota. In adjusted models, MDI and PDI scores were not associated with alpha diversity indices. However, in adjusted models, MDI and PDI scores were positively associated with the first coabundance factor, which captured positive loadings for the genera Faecalibacterium, Sutterella, and Clostridium cluster XIVa. For an interquartile range increase in the first coabundance factor, MDI scores increased by 3.9 points [95% confidence interval (CI): 0, 7.7], while PDI scores increased by 8.6 points (95% CI 3.1, 14). Our results highlight the potential for gut microbial compositional characteristics to be important correlates of children's Bayley Scales performance at 36 months of age.

Correlation between children's gut microbiota and maternal/child characteristics. For coabundance factor 1, children born by Cesarean birth (C-section) had lower values, compared to children born vaginally, while no differences were observed in alpha diversity measures (Supplementary Tables S1, S2). Neonatal gut microbiota differ between children born vaginally and by C-section 16 ; however, the long-term impacts vary 17,18 . This may be due to other factors that impact children's later gut microbiota composition (e.g., infant feeding practices and child's diet) 5 .
For coabundance factor 2, children consuming fish within the previous 24 h had higher values, compared to children who did not (Supplementary Table S1). Coabundance factor 2 captured positive loadings for the genera Roseburia, Coprococcus, and Blautia (Table 2). Fish tissue is enriched in omega-3 fatty acids, including docosahexaenoic acid (DHA) 19 , which has been associated with enrichment of these taxa or taxa from the same family (Lachnospiraceae) 20,21 .
The second and third coabundance factors and/or alpha diversity indices were lower among children who had illnesses in the previous 12 months, including fever and respiratory illness, compared to children who did not (Supplementary Tables S1, S2), which was consistent with prior studies concerning childhood illnesses and reduced gut microbiota diversity 4,5 . The third co-abundance factor was also positively correlated with both gestational age and higher household income (Supplementary Tables S1 and S3), although the latter two variables were not correlated to each other (Wilcoxon rank sum test, p = 0.73).
Multivariable regression models using coabundance factors. In both MDI and PDI models, the first coabundance factor was positively correlated with both outcome measures (Table 3). In adjusted models, for an interquartile range increase in the first coabundance factor 1 (loading: 2251), MDI scores increased on average by 3.9 points [95% Confidence Interval (CI): 0, 7.7], while PDI scores increased on average by 8.6 points (95% CI 3. 1,14). Compared to other covariates, the first coabundance factor explained the largest proportion of the variance in PDI scores (24%), while in the MDI model, the first coabundance factor explained 12% of the variance. In contrast, in the MDI and PDI models, the second coabundance factor explained < 1% and 9% of the variance, respectively, while the third coabundance factor explained 7% and 1% of the variance, respectively. www.nature.com/scientificreports/ www.nature.com/scientificreports/ Child's weight-for-age (z-score) was associated with higher MDI scores (β = 3.2 points, 95% CI: 0.76, 5.7), while each 1-month increase in breastfeeding was associated with a 0.8-point (95% CI: − 0.01, 1.6) increase in the MDI score (Table 3). In contrast, the strongest predictors of PDI scores were child sex, with boys having lower scores compared to girls (β = − 11 points, 95% CI: − 18, − 3.1), while child consumption of fish in the previous 24 h was associated with higher PDI scores (β = 12 points, 95% CI: 2.2, 21) (Table 3). Table S3), and in adjusted regression models (Supplementary Table S4).

Multivariable regression models using alpha diversity measures and Maaslin2. Alpha diversity measures were not strongly correlated with Bayley MDI and PDI scores in bivariate analyses (Supplementary
Using Maaslin2, in adjusted models, several taxa were correlated with Bayley MDI and PDI scores, although q-values were > 0. 25 (Supplementary Table S5). The genus Faecalibacterium was positively correlated with both outcome measures, while the genus Gemmiger was positively correlated with PDI scores. Both taxa were among the dominant genera in the first coabundance factor (positive loadings: 0.59 and 0.58, respectively) ( Table 2), which was also positively correlated with both MDI and PDI scored in adjusted models ( Table 3). The genus Lachnospiraceae incertae sedis was negatively correlated with MDI scores, while the genus Flavonifractor was positively correlated with MDI scores (Supplementary Table S5). For both genera, the same directions of association were observed with coabundance factor 1 (negative and positive, respectively), although they were not among the dominant genera (loadings: -0.18 and 0.09, respectively) (  Table S7). For both analyses, the same patterns of association between Bayley scores and coabundance factors were observed as in our main analysis (Table 3).

Discussion
At 36 months of age, children's gut microbiota composition was associated with Bayley scores. In adjusted models, the first bacterial coabundance factor (dominated by positive loadings for the genera Sutterella, Oscillibacter, and Alistipes, and negative loadings for the genera Anaerostipes, Clostridium cluster XVIII, and Blautia) was positively correlated with both MDI and PDI scores (Table 3), and explained a relatively large proportion of the variability in these scores in adjusted models (12% and 24%, respectively). Our results highlight the potential for microbial compositional characteristics to be important correlates of Bayley Scales performance in our small cross sectional assessment.
Findings in our study differed from previous studies investigating children's gut microbiota composition and neurodevelopment [6][7][8] . In North Carolina, USA, children's gut microbiota were clustered into three groups; the cluster represented by Faecalibacterium had the lowest cognitive scores (assessed using the Mullen Scales of Early Learning), while measures of alpha diversity (observed number of species, Chao1, Faith's phylogenetic diversity, and the Shannon index) were also associated with lower cognitive scores 7 . In our study, Faecalibacterium was associated with higher MDI and PDI scores ( Table 3, Supplementary Table S4), while alpha diversity indices were not strongly correlated with MDI and PDI scores (Supplementary Table S3). In Massachusetts, Missouri, and California, USA, the bacterial coabundance factor with positive loadings for Lachnospiraceae was associated with poorer communication and personal and social skills in 3-year-olds (assessed using the Ages and Stages Questionnaire, 3rd edition) 6 . This differed from our study, in which the first coabundance factor contained both positive and negative loadings for genera of the Lachnospiraceae family, and this coabundance factor was correlated with better MDI and PDI scores (Table 3). Likewise, using Maaslin2, members from the Lachnospiraceae family (Lachnospiraceae incertae sedis and Clostridium cluster XlVb) were associated with both lower and higher MDI scores, respectively (Supplementary Table S5). In Ohio, USA, alpha diversity (phylogenetic diversity and the Shannon index) was associated with children's temperament, which was assessed by maternal report using the Early Childhood Behavior Questionnaire 8 . As noted above, Bayley Scales and alpha diversity indices were not strongly correlated in our study (Supplementary Tables S3, S4). Differences between our findings and results from other studies may reflect temporal development of the children's microbiota 5 , as well as potential differential sensitivity to the microbiome during different periods of neurodevelopment. In prior studies, children's gut microbiota were assessed at 6 months 6 , 1 year 7 , and 2 years 8 , while in our cross-sectional study, children's gut microbiota were assessed at 36-37 months, which is considered a stable phase in microbiome progression 5 .
One mechanism by which the gut microbiota is hypothesized to influence the central nervous system, is through the production of metabolites, including short-chain fatty acids (SCFAs) 22,23 . SCFAs enhance the integrity of the intestinal epithelial layer 24 , and also travel to other organs, including the brain 22,25 . In germ-free mice, the lack of a normal gut microbiota was associated with increased permeability of the blood-brain barrier 22 . Moreover, colonization of germ free mice by bacterial strains that produced SCFAs, including butyrate or acetate/ propionate, reversed this trend and decreased permeability of the blood-brain barrier 22 .
Coabundance factor 1 captured positive loadings for the genus Faecalibacterium (loading: 0.59), which is considered one of the dominant butyrate producers and has been linked to anti-inflammatory effects 26 . Coabundance factor 1 also captured positive loadings for the genus Clostridium cluster XIVa (loading: 0.59), which includes a large number of butyrate-producing bacteria 27,28 , and positive loadings for the genus Alistipes (loading: 0.79), a propionate-producer 29 . Blautia spp. produce propionate 29 ; however, this genus was negatively correlated with coabundance factor 1 (loading: − 0.66). Using Maaslin2, several butyrate-and propionate-producing bacteria were positively correlated with MDI and/or PDI scores, including Faecalibacterium 26  www.nature.com/scientificreports/ These results suggest the production of the SCFAs, including butyrate and/or propionate, may contribute positively to children's Bayley scores in our cohort. However, as noted above, this was not consistent across all components of coabundance factor 1. Furthermore, coabundance factor 2 captured positive loadings for Roseburia, Coprococcus, and Blautia (range for positive loadings: 0.45, 0.82) ( Table 2), which are also known to produce butyrate and/or propionate 29 , and coabundance factor 3 captured positive loadings for the genera Roseburia and Butyricicoccus (loadings: 0.52, 0.64), which are both known to produce butyrate 29 . But coabundance factors 2 or 3 were not associated with Bayley scores (Table 3 and Supplementary Table S3), suggesting the contributions of butyrate and/or propionate production to associations were uncertain, and other factors likely contributed.
The highest positive loadings for coabundance factor 1 were for Alistipes, Sutterella, and Oscillibacter (positive loadings: 0.79, 0.85, and 0.88, respectively). In other studies, the same genera were adversely associated with cognitive outcomes [32][33][34][35] , which differed from the direction of association observed in our study. Sutterella spp. are considered a common intestinal commensal among adults and children, which decrease in abundance between the small intestine and the colon [36][37][38] . Studies concerning the role of Sutterella among children have focused on autism spectrum disorder (ASD). Compared to controls, among children with ASD, Sutterella spp. were higher in fecal samples 34 , and in ileal and cecal biopsies 35 . Similarly, clades within the genera Oscillibacter and Alistipes were enriched among adults with depression, compared to healthy controls 32,33 . Oscillibacter may contribute to depression due to production of valeric acid, which resembles the neurotransmitter GABA (gammaaminobutyric acid) 33 . Disruption in GABA signaling has been linked to neurodevelopmental disorders 39 . Thus, in previous studies, all three genera were associated with poorer neurobehavioral function (childhood ASD and adult depression) [32][33][34][35] , whereas, in our analysis, they were associated with better neurocognitive function, i.e., higher MDI and PDI scores in association with increasing values for coabundance factor 1 ( Table 3). This suggests the roles of these taxa may potentially vary by population characteristics (e.g., age, neurologic status) or by timing of microbiota assessment, and have yet to be elucidated. Lastly, in a cross-sectional analysis, there is the potential for reverse causation which may lead to divergent findings across studies, depending on how (or if) the outcome of interest impacts microbiota characteristics.
Although our study has many strengths, there are some limitations to note. Our sample size was limited; however, this study was exploratory. Despite the small sample size, our results suggest that gut ecology is an important correlate of children's neurodevelopment in our cohort. There were some differences between participants who collected their child's stool sample and those that did not (Table 1), which could have contributed to selection bias. The possibility of residual or unmeasured confounding remains, although we adjusted for many potential confounders. For example, we did not have information on potentially relevant factors, such as antibiotic treatment, which is known to impact the gut microbiota 4 . This was a cross-sectional study, and therefore it was not possible to determine the direction of association between children's gut microbiota and outcome measures. In addition, child's stool collection (at 3 years of age) occurred when microbiome development is considered more stable 5 . One potential ramification of measuring the microbiome when it is relatively stable is the underestimation of previous variation in the microbiome that could have impacted neurodevelopment. This is a form of measurement error and, as it is likely to be non-differential error, it would bias findings to the null. Thus, this timing of gut microbial assessment may lead to an underestimation of the "true" association. Longitudinal studies are needed to investigate whether there is a developmental window for the gut microbiome, which is most predictive for children's neurodevelopment.
In conclusion, in our cross sectional analysis of 36-month old children living in rural China, children's gut microbiota explained a large portion of variability in Bayley Scales (MDI and PDI scores). Results from this study suggest that gut microbiota may play an important role in children's neurodevelopment among typically developing children.

Methods
Recruitment. The study participants were residents of Daxin County, Guangxi province, China, which is predominantly rural, with a population of 359,800 (http:// www. daxin. gov. cn/). Approximately 50,000 residents live in the small rural town of Daxin, and most of the remaining are rice farmers and live outside of town.
Between May 2013 and March 2014, 408 mother-infant pairs were recruited at the Maternal and Child Health Hospital in Daxin County. Eligible mothers were in good general health, resided in Daxin County during the previous 3 months, and planned to remain in Daxin County for the next 12 months. Mothers were enrolled up to 4 weeks before delivery or within one week postpartum. Of 408 mothers enrolled in the study, 10 mothers were subsequently excluded because they actually lived outside Daxin County (n = 3), gave birth to twins (n = 1), or data collection was incomplete (n = 6), resulting in a cohort of 398 mother-infant pairs. Participants returned between May 2016 and March 2017 for their child's 36-month assessment.
The research protocol was reviewed and approved by the Institutional Review Boards at Oregon State University (USA), the University of South Carolina (USA), and Xinhua Hospital (China). Mothers (or the child's primary caregiver) provided written informed consent prior to their (and their child's) participation. All methods were carried out in accordance with relevant guidelines and regulations.
Questionnaires. After enrollment, while in the hospital, mothers filled out a detailed questionnaire ascertaining demographics, education, occupation, household income, pre-pregnancy weight and height (for body mass index calculation), delivery mode, and their child's sex. Gestational age was abstracted from the medical record based on maternal reported date of the last menstrual period and date of delivery. At 36 months, parents or caregivers completed a questionnaire regarding their child's health, breastfeeding duration, preschool attendance, primary caregiver, and their child's food intake during the previous 24 h. To characterize the child's health, caregivers were also asked whether the child had any of the following seven illnesses or symptoms in the previ- www.nature.com/scientificreports/ ous 12 months: upper respiratory illness, lower respiratory illness, difficulty breathing, diarrhea, vomiting, fever, and rash; the frequency with which the illness or symptom occurred; and whether the child saw a doctor for the given illness or symptom. Children's weight and height were measured by hospital staff using a digital standing scale/stadiometer (Model # HCS-100-RT, Jiangsu, China). Z-scores for weight-for-age, height-for-age, and weight-for-height were calculated using the 2006 World Health Organization (WHO) child growth standards 40 , with the R package "anthro" (21 May 2020). Weight-for-age (z-score) reflects short-term nutritional deficiencies, while long-term deficiencies reduce stature, which are reflected in height-for-age and weight-for-height z-scores 41 .
Neurodevelopmental assessment. At 36 months, children's neurodevelopment was assessed using the BSID-II, which was translated into Chinese. The BSID-II yields two main indices of performance: the Mental Developmental Index (MDI) and the Psychomotor Developmental Index (PDI) 12 . These indices are age-standardized to a mean of 100 and standard deviation of 15, based on an English language-speaking U.S. reference population 12 . The BSID-II underwent extensive pre-testing, and all children were evaluated by one doctor (coauthor YN), who spoke the local dialect. The evaluator was trained on-site in Bayley administration by co-author EPT, a developmental psychologist. Examiner reliability was assessed throughout the follow-up period by videotaping five children (of 196 children = 2.6%). Videotapes for both MDI and PDI sections were evaluated and rescored by co-authors EPT and FJB (also a developmental psychologist), and differences in scoring were minor.
Analysis of children's gut microbiota. Participants collected their child's stool sample at home using a sterile cotton swab (Becton Dickinson Cat #220135, Franklin Lakes, NJ). The cotton swab was removed from the tube, and then used to swipe toilet paper containing the child's feces.  42 . The concentration of extracted DNA was determined using a BioDrop (Biochrom Ltd., Cambridge, UK); its integrity and size was checked by 0.8% agarose gel electrophoresis containing 0.5 mg/ml ethidium bromide. Hypervariable region V3-V4 amplicons from the 16S rRNA gene were sequenced by Illumina MiSeq 2 × 300 bp paired-end sequencing (https:// suppo rt. illum ina. com/ docum ents/ docum entat ion/ chemi stry_ docum entat ion/ 16s/ 16s-metag enomic-libra ry-prep-guide-15044 223-b. pdf). Detailed methods are in the Supplementary Information.
Quality-filtered reads were delineated into unique sequences and then sorted by decreasing abundance, and singletons were discarded. Using previously described methods 43 , operational taxa units (OTUs) were clustered de novo with Uparse at 97% similarity level 44 . Reference-based chimera detection was performed using UCHIME 45 against the RDP (Ribosomal Database Project) classifier training database (v9) 46 . The OTU table was finalized by mapping quality-filtered reads to the remaining OTUs with the Usearch 47 global alignment algorithm at a 97% cutoff. Sequence data were rarefied to 13,000 reads per sample (1000 permutations) to avoid bias caused by the difference in sequencing depth. Alpha diversity (within person diversity) was assessed using the number of observed OTUs, Shannon's diversity index, Faith's Phylogenetic Diversity, and Pielou's measure of evenness.
Coabundance groupings. Coabundance groupings (or factor representations) of bacterial taxa within the children's gut were computed, similar to previously described methods 6 . Briefly, Spearman's rank correlation coefficients were determined for the 25 most abundant taxa (based on the median relative abundance), using the Compositionality Corrected by Permutation and Renormalization ("ccrepe") package in R 48 . Spearman's correlation matrix was used for principal components analysis (with varimax rotation), using the "psych" package in R 49 . This analysis revealed a 3-factor solution based on a scree plot. Factor loadings were applied to 16S rRNA sequencing counts to compute an individual's factor score, for each bacterial coabundance group. Data analysis. Differences in maternal/child characteristics between participants who provided a stool sample and those who did not were investigated using a Chi-squared test, Fisher's exact test, or Wilcoxon rank sum test. Bivariate associations between Bayley scores, children's gut microbiota (coabundance scores and alpha diversity measures), and maternal/child characteristics were then investigated using a Spearman's correlation, Wilcoxon rank sum test, or Kruskal-Wallis test. For bivariate associations and multivariable regression models, the number of levels for some categorical variables were combined: mothers completing high school and university were combined into one level (maternal education: < high school, ≥ high school), mothers who were workers (civil servant, white-collar worker, skilled worker, unskilled worker, and shopkeeper) or listed "other" as the type of work were combined into one category for non-farmers (maternal occupation: farmer, non-farmer), households earning 2000-5000 RMB per month and ≥ 5000 RMB per month were combined into one level (household monthly income: < 2000 RMB, ≥ 2000 RMB), and children cared for by their mothers or fathers were combined into one level, and were compared to children cared for by a grandparent (child's primary caregiver: mother/ father, or grandparent).
Multivariable linear regression models were constructed to examine the relationships between Bayley scores (MDI and PDI) and children's gut microbiota. Primary analyses included all three coabundance factors as independent predictors. Table 1 includes all potential covariates that were available and considered for models based on evidence that they were independently associated with children's neurodevelopment and/or the gut microbiota www.nature.com/scientificreports/ in previous studies [4][5][6][7][8][9][10][11][15][16][17][18][19][20][21] . From this list, we developed a core model based on priors that included child's sex, child's age at BSID-II assessment, breastfeeding duration, and maternal age and education. The subset of additional Table 1 covariates that were included in the final models were chosen based on at least one of the following: their relationship with outcome measures; their relationship with coabundance factors; added variable plots; evidence of confounding of effect estimates for coabundance factors; and comparison of the Akaike Information Criterion between models with/without a given covariate 50 . From Table 1, child's rice consumption was the only variable not considered for inclusion because all children ingested rice or rice porridge in the previous 24 h. Associations between Bayley scores and gut microbiota may differ between individual taxa and coabundance factors. Therefore, secondary analyses investigated the associations between covariate-adjusted Bayley scores (modeled as residuals) and individual taxa (genus-level), using Microbiome Multivariable Association with Linear Models (Maaslin2) in R 51 . For Maaslin2, individual taxa were included if measured in at least 20% of children's gut microbiota. Bacterial taxa were transformed using the default transformation in Maaslin2 (arcsinesquare-root transformation), which has been applied to proportional data 52 . Secondary analyses also investigated the associations between Bayley scores and each alpha diversity measure in multivariable regression models.
Among the participants included in this analysis, just one variable (household income) had missing values, and missing data were imputed using multiple imputation based on the multivariate normal distribution 53 , as previously described 10 . Associations were checked with and without imputed data, and no differences were observed. As a sensitivity analysis, regression models using coabundance factors were re-evaluated using the Bayley raw scores in place of the age-standardized scores. One child was born pre-term (35.3 weeks), and as a sensitivity analysis, models using coabundance factors were re-evaluated including only children with term births (≥ 37 weeks). Assumptions for model residuals were checked (no evidence of non-linearity, constant variance, normal distribution) using appropriate plots. Potential influential observations were checked using Cook's distance; however none were observed.
For bivariate analyses and regression models, an alpha level of 0.05 was used as a guide for significance. To minimize the false discovery rate (FDR), we used dimensionality reduction (bacterial coabundance groups) in our primary analyses to characterize the gut microbiota, and thus q-values were not calculated. Using Maaslin2, q-values (FDR-corrected p-values) were calculated 46 , and q ≤ 0.25 was used as a guide for significance. Statistical analyses were performed using Stata (Version 9.2, College Station, TX, USA), and the R-platform (Version 4.0.2, 06 June 2020).