An overview of the methods is presented in Fig. 1. The RESONANCE cohort and NHBCS are part of ECHO, a consortium that harmonizes data across cohorts to probe scientific questions regarding children’s environmental health that are unanswerable by a single cohort [27, 28]. Background on each of the cohorts, including eligibility criteria, is provided in the Supplemental Methods.

Diagram of Study Methods. Fecal samples were collected from RESONANCE participants between six weeks and two years postpartum and from New Hampshire Birth Cohort Study (NHBCS) participants at six weeks and/or one year. In both cohorts, DNA was extracted from stools and underwent shotgun metagenomic sequencing. In both cohorts, assessment of autism-spectrum disorder behaviors with the Social Responsiveness Scale, 2nd edition, was completed after age 3 years. Additionally, a subset of NHBCS fecal samples underwent nuclear magnetic resonance spectroscopy to quantify small molecules. These data were integrated with sociodemographic data into statistical models

Fecal samples were collected in tubes containing an RNA stabilizer from a subset of child participants for microbiome sequencing and analysis. In the NHBCS, fecal samples (n = 389) were collected from 212 children at approximately 6 weeks and 12 months of age, stored in a polyethylene bag, and frozen in home freezers or directly transported on ice to maternal postnatal follow-up visits. Within 24 h of receipt, the samples were aliquoted and frozen at 80 °C. We restricted this analysis to samples collected before four years of age (n = 92). DNA was extracted from the stool using the Zymo DNA extraction kit (Zymo Research, Irvine, California) and prepared for metagenomic sequencing using Illumina NextSeq at Marine Biological Laboratory in Woods Hole, Massachusetts. In the RESONANCE cohort, stool samples were collected by parents in OMR-200 tubes (OMNIgene GUT, DNA Genotek, Ottawa, Ontario, Canada), immediately stored on ice, and brought within 24 h to the lab in Rhode Island where they were immediately frozen at − 80 °C. Stool samples were not collected if the subject had taken antibiotics within the last 2 weeks. DNA extraction was performed at Wellesley College (Wellesley, Massachusetts). Nucleic acids were extracted from stool samples using the Rneasy PowerMicrobiome kit automated on the QIAcube (Qiagen, Germantown, Maryland), excluding the DNA degradation steps. Shotgun metagenomic sequencing was performed on extracted DNA at the Integrated Microbiome Resource (IMR, Dalhousie University, Nova Scotia, Canada).

For the NHBCS, sequencing libraries were prepared using Nugen’s Ovation Ultralow V2 protocol. The extracted DNA samples were sheared to a mean insert size of 400 bp using a Covaris S220 focused ultrasonicator, resulting in an average of 30 million reads per sample. For the RESONANCE cohort, a pooled library (max 96 samples per run) was prepared using the Illumina Nextera Flex Kit for MiSeq and NextSeq from 1 ng of each sample. Samples were then pooled onto a plate and sequenced on the Illumina NextSeq 550 platform using 150 + 150 bp paired-end “high-output” chemistry, generating ~ 400 million raw reads and ~ 120 Gb of sequence per plate. For both cohorts, the BioBakery v3 pipeline was used to process the raw sequences [29]. KneadData (v0.7.7) was used to trim adapters from raw sequence reads and remove reads matching a human genome reference. Metagenomic Phylogenetic Analysis (MetaPhlAn v3.0.7) was used to generate taxonomic profiles by aligning reads to a reference database of marker genes (mpa_v30_CHOCOPhlAn_201901) [29]. Finally, the Human Microbiome Project Unified Metabolic Analysis Network (HUMAnN v3.0.0a4 with databases at v201901b) was used to functionally profile metagenomes [29] .

In the NHBCS, stool samples were collected for metabolomics analysis in trace element free tubes and sent to the Eastern Regional Comprehensive Metabolomics Research Core at the Research Triangle Institute for nuclear magnetic resonance (NMR) analysis. NMR spectra of fecal samples were acquired on a Bruker 700 MHz NMR spectrometer (Billerica, Massachusetts). The relative concentrations of library-matched metabolites associated with host and gut microbe co-metabolism [30,31,32] were obtained using the Chenomx NMR Suite 8.1 Professional (Chenomx Inc., Edmonton, Alberta, Canada) [33, 34]. Relative concentrations were log2-transformed for analysis.

Caregivers assessed the participant’s social skills using the Social Responsiveness Scale, 2nd edition (SRS-2), which captures ASD-related social behaviors in multiple domains, in both cohorts [35]. Standardized total T-scores have a mean (standard deviation, SD) of 50 (10), with higher scores indicative of social behaviors associated with ASD. In the NHBCS, caregivers completed the preschool version of the SRS-2 when participants were ~ 3 years of age and the school-age version when the participants were 5 years of age. If a preschool assessment was not available, the participant was included with their school-age score. In the RESONANCE cohort, caregivers completed the SRS-2 preschool version for participants 2.5–4 years of age and the school-age version for participants 4–19 years of age.

Covariates were selected a priori based on their potential to confound the association between the microbiome and SRS-2 scores or as predictors of the outcome. At recruitment, study staff implemented questionnaires to collect data on maternal smoking during pregnancy (any vs. none), parity (parous vs. nulliparous), parental age (continuous), and sociodemographic characteristics, including maternal education (any higher education vs. none) and marital status. Characteristics of the birth (delivery mode: C-section vs. vaginal, peripartum antibiotic exposure: any vs. none, gestational age at delivery: continuous, and child’s sex) were abstracted from medical records (NHBCS) or questionnaires (RESONANCE). Repeated questionnaires throughout childhood captured duration of breastfeeding (continuous) and breastfeeding status at the time of stool sample collection (exclusive vs. any formula). Exact age at SRS-2 implementation was calculated as the difference between the date of birth and the date of survey completion (continuous). To reduce selection bias and improve power, we imputed missing covariates using the mice R package [36,37,38]. Details are available in the Supplemental Methods.

Because many NHBCS participants provided two stool samples, we analyzed early (< 6 month) and late (6 month–2 year) samples separately. Additionally, we examined NHBCS and RESONANCE samples separately. We tested the overall association between bacterial species or functions and SRS-2 scores with adonis2, a PERMANOVA method that allows for quantification of variable influence (R2) and significance (marginal p-value) in community structure (Bray–Curtis distances) [39,40,41]. A p-value less than 0.05 was considered significant. Bacterial diversity was quantified with Shannon and Inverse Simpson indices. We linearly regressed diversity against SRS-2 scores, adjusting for the previously mentioned covariates, considering a p-value less than 0.05 significant. We tested whether models built on data from the NHBCS could predict SRS-2 scores in the RESONANCE cohort using the caret package in R [42]. Briefly, we used the predict function in the stats R package to predict SRS-2 scores for RESONANCE observations based on their microbiome and covariates and the model constructed using data from the NHBCS. Functions from the caret package were then used to calculate correlations and errors between actual and predicted scores. To screen for bacterial species associated with SRS-2 scores, we used microbiome multivariable association with linear models (MaAsLin2), which models bacterial relative abundance as the outcome [43]. Associations with a false discovery rate (FDR) of q < 0.1 were modeled in linear regression with bacterial relative abundance predicting SRS-2 scores, adjusting for covariates. Given the sex-ratio differential of ASD presentation, we explored whether relationships between the gut microbiome and behavior were different in biological males and females with interaction models [44].

Potentially neuroactive microbial gene sets were acquired from Supplementary Dataset 1 from Valles-Colomer et al. [45] Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthologs (KOs) from this dataset were mapped to UniRef90s using the mapping files included with HUMAnN3. Pearson correlations for each measure (SRS-2 scores or residuals from a complete linear model regressing SRS-2 scores against all covariates) and each UniRef90 gene family found in ≥ 4 subjects were calculated. We performed Mann–Whitney U tests (implemented in the HypothesisTests.jl package) on the correlations of each gene family in a neuroactive gene set against the correlations of all gene families not found in the gene set [46]. FDR q < 0.1 was the threshold for significance [47].

To test if the capacity of the gut microbiome to have neuroactive genes associated with SRS-2 scores related to the metabolic profiles of infants, we regressed SRS-2 scores against relative concentrations of fecal metabolites identified as significant from the Mann–Whitney U test that were also annotated in our dataset. Contribution of specific taxa to the relative abundance of neuroactive gene pathways was determined using stratified MetaPhlAn tables. All models were run with R (4.1.1), RStudio (v1.4.1717), and Julia (v1.6.1). Code and package versions are available at https://github.com/HEL548/ECHOSRS2.