Autism-Related Dietary Preferences Mediate Autism-Gut Microbiome Associations
Published:November 11, 2021DOI:https://doi.org/10.1016/j.cell.2021.10.015
Highlights
Summary
There is increasing interest in the potential contribution of the gut microbiome to autism spectrum disorder (ASD). However, previous studies have been underpowered and have not been designed to address potential confounding factors in a comprehensive way. We performed a large autism stool metagenomics study (n = 247) based on participants from the Australian Autism Biobank and the Queensland Twin Adolescent Brain project. We found negligible direct associations between ASD diagnosis and the gut microbiome. Instead, our data support a model whereby ASD-related restricted interests are associated with less-diverse diet, and in turn reduced microbial taxonomic diversity and looser stool consistency. In contrast to ASD diagnosis, our dataset was well powered to detect microbiome associations with traits such as age, dietary intake, and stool consistency. Overall, microbiome differences in ASD may reflect dietary preferences that relate to diagnostic features, and we caution against claims that the microbiome has a driving role in ASD.
Graphical abstract
01231-9&id=fx1.jpg){kind=link}
Introduction
There is considerable interest in the relationship between the gut microbiome and autism spectrum disorder (ASD)—a neurodevelopmental condition characterized by social and communication difficulties and restricted and repetitive behaviors as well as unusual sensory responsiveness. This interest is driven by high rates of co-occurring gastrointestinal symptoms (Chaidez et al., 2014; Kohane et al., 2012; McElhanon et al., 2014; Vargason et al., 2019) and by animal studies suggesting that the microbiome causally contributes to ASD-related behavioral traits (Buffington et al., 2016; Hsiao et al., 2013; Sharon et al., 2019). There is also increasing evidence implicating a brain-gut-microbiome axis in a range of gastrointestinal (e.g., irritable bowel syndrome, inflammatory bowel disease) and neuropsychiatric conditions (e.g., depression) (Valles-Colomer et al., 2019).Despite this intense interest and the recent progression of fecal microbiota transplantation studies in ASD to phase I clinical trials (Kang et al., 2019; Kang et al., 2017), the contribution of the microbiome to ASD is unconvincing. Existing human ASD microbiome studies have yielded inconsistent results, likely due to underpowered study designs (e.g., n < 50) (De Angelis et al., 2013; Kang et al., 2018) that are prone to bias (Rothschild et al., 2020), variable use of appropriate multiple-testing correction (De Angelis et al., 2013; Parracho et al., 2005; Williams et al., 2011), inconsistent use of compositionally aware statistical analysis (De Angelis et al., 2013; Finegold et al., 2010; Williams et al., 2011), and large site effects that may represent non-uniform processing (Fouquier et al., 2021). Indeed, many of the larger studies with greater than 100 participants have found no ASD-microbiome associations (Gondalia et al., 2012; Son et al., 2015). Furthermore, with few exceptions (Berding and Donovan, 2018; Son et al., 2015), most ASD microbiome studies have not comprehensively considered important microbiome confounders (Bokulich et al., 2016; Falony et al., 2016; Rothschild et al., 2018; Vandeputte et al., 2016) (via design or statistical analysis) including sex (Adams et al., 2011; Wang et al., 2011), age (Adams et al., 2011), diet (Adams et al., 2011; De Angelis et al., 2013; Finegold et al., 2010; Strati et al., 2017; Wang et al., 2011), stool consistency as measured using the Bristol Stool Chart (rBSC), or co-occurring gastrointestinal complaints and medications (Wang et al., 2011). Notably, ASD has a male bias, and children on the autism spectrum often have less-diverse diets (Panossian et al., 2020; Schreck and Williams, 2006) (due to selective eating and allergy concerns) and co-occurring gastrointestinal conditions (Chaidez et al., 2014; Kohane et al., 2012; McElhanon et al., 2014) and are more likely to be prescribed antibiotics (Niehus and Lord, 2006). Thus, it is unclear whether existing reports of microbiome associations in ASD may be explained by confounding variables. There are also technical limitations: most studies have used 16S rRNA sequencing (Ho et al., 2020) (with the exception of a few [Dan et al., 2020; Wan et al., 2021; Wang et al., 2019]), which provides lower taxonomic resolution and limited functional information about the microbiome.Within the limitations of these studies, a recent meta-analysis reported a small number of taxa with consistent evidence of association with ASD diagnostic status; namely, the genus Prevotella, the phylum Firmicutes, clusters of the order Clostridiales, and species of Bifidobacterium (Ho et al., 2020). Nonetheless, to our knowledge, the overall relationship between the microbiome and ASD diagnosis has not been quantified. In fact, it appears that interest in the ASD microbiome outstrips what the primary evidence base warrants: in April 2021, a PubMed search ([autism [Title]] AND [microbiome[Title]] OR [microbiota[Title]]) identified 56 review articles, compared to only 26 primary research studies, included in the latest meta-analysis (Ho et al., 2020). More broadly, a meta-analysis of human microbiota-associated animal studies has raised concerns that the sheer extent of positive findings is implausible, suggesting that causal inference has been premature, over-stated, and suffers from publication bias (Walter et al., 2020).Here, we present a stool metagenomics study of 247 children from the Australian Autism Biobank (AAB) (Alvares et al., 2018) and the Queensland Twin Adolescent Brain (QTAB) project, for which extensive phenotype data were available (Figures 1A and 1B; Table 1): demographics, dietary data from the Australian Eating Survey (from which we derived dietary diversity and principal components from percentage energy measures, hereafter referred to as dietary PC1-3) (see STAR Methods; Figure 1B), stool consistency, detailed psychometric testing, and genome-wide single nucleotide polymorphism (SNP) genotypes. We find that the stool microbiome captures negligible variation in ASD diagnosis, whereas there were large and significant associations with variance in dietary traits, stool consistency, and age. Instead, our results suggest an alternative model whereby genetic and phenotypic measures of the autism spectrum (including restricted and repetitive behaviors, social affect, and higher sensory sensitivity) promote a less diverse diet, which reduces microbiome diversity and is associated with looser stool consistency.
01231-9&id=gr1.jpg){kind=link}
Table 1Demographic summary of participants in the metagenomics study
| Phenotype | AAB+QTAB microbiome cohort | UNR sub-groups | Entire AAB | |||||
|---|---|---|---|---|---|---|---|---|
| ASD | SIB | UNR | UNR (AAB) | UNR (QTAB) | ASD | SIB | UNR | |
| N | 99 | 51 | 97a | 48 | 49 | 1168 | 262 | 149 |
| Age | 8.7 (3.8) | 8.0 (4.3) | 9.2 (3.7) | 6.3 (3.0) | 12.1 (1.0) | 7.7 (3.9) | 8.2 (4.2) | 6.5 (3.4) |
| Male % | 73% | 57% | 54% | 50% | 59% | 60% | 40% | 39% |
| BSC | 3.56 (1.31) | 3.84 (1.25) | 3.37 (0.96) | 3.35 (1.08) | 3.39 (0.83) | 3.51 (1.02) | 3.56 (0.97) | 3.46 (0.92) |
| rBSC | 3.48 (0.96) | 3.71 (0.94) | 3.40 (0.84) | 3.39 (0.91) | 3.41 (0.78) | 3.58 (1.42) | 3.64 (1.23) | 3.44 (1.18) |
| dietary PC1 | −0.15 (1.85) | 0.13 (1.33) | 0.09 (0.93) | 0.07 (1.01) | 0.10 (0.87) | – | – | – |
| dietary PC2 | −0.20 (1.60) | −0.14 (1.36) | 0.28 (1.10) | 0.23 (1.06) | 0.32 (1.15) | – | – | – |
| dietary PC3 | −0.30 (1.73) | 0.28 (0.67) | 0.16 (0.66) | −0.11 (0.61) | 0.43 (0.6) | – | – | – |
| ARFS | 28.6 (11.5) | 34.0 (9.8) | 33.8 (9.7) | 32.9 (9.2) | 34.6 (10.2) | – | – | – |
| Taxonomic Shannon index | 3.63 (0.44) | 3.59 (0.45) | 3.76 (0.38) | 3.59 (0.40) | 3.92 (0.28) | – | – | – |
| Dietary Shannon index | 4.61 (0.06) | 4.63 (0.05) | 4.65 (0.04) | 4.65 (0.04) | 4.66 (0.03) | – | – | – |
| Genome-wide SNP data N | 91 | 43 | 84 | 48 | 49 | 887 | 218 | 116 |
| ASD PGS | 0.18 (1.10) | 0.22 (0.95) | 0.028 (0.72) | 0.06 (1.00) | −1e-6 (1.9e-7) | 0.060 (1.00) | −0.085 (0.96) | 0.033 (1.07) |
| ASD-ADHD-TS PGS | 0.10 (1.10) | 0.17 (0.82) | 0.04 (0.61) | 0.09 (0.87) | −1.1e-6 (1.6e-7) | 0.047 (0.98) | −0.079 (0.90) | 0.012 (0.88) |
| Neuroticism PGS | −0.038 (1.00) | −0.086 (0.94) | 0.041 (0.68) | 0.084 (0.98) | −4.1e-5 (4.2e-5) | 0.040 (1.00) | 0.025 (0.95) | 0.145 (0.93) |
| IQ-DQ composite score | 85 (24) | 100 (17) | 100 (15) | 100 (15) | 100 (15) | 79 (24) | 101 (13) | 103 (13) |
| CSHQ raw score | 43 (11) | 37 (6.1) | 38 (7.1) | 38 (7.1) | – | 44 (10) | 39 (9) | 37 (8) |
| ADOS2/G comparison score | 6.7 (2.1) | – | – | – | – | 6.7 (2.0) | – | – |
| ADOS2/G RRB CSS score | 7.0 (2.3) | – | – | – | – | 6.8 (2.2) | – | – |
| ADOS2/G social affect CSS score | 6.8 (2.0) | – | – | – | – | 6.8 (2.0) | – | – |
| SSP Sensory raw score | 41 (12) | – | – | – | – | 33 (8) | – | – |
| SRS T-score | 74 (9.3) | 46 (9.1) | 51 (9.3) | 51 (9.3) | – | 78 (10) | 50 (12) | 50 (9) |
First 3 columns show the comparison across ASD/SIB/UNR groups used in this analysis. Columns 4 and 5 show breakdown of UNR group into AAB and QTAB cohorts. Columns 6–8 show characteristics of the broader AAB cohort. BSC: Bristol Stool Chart (scale 1–7), rBSC: regrouped BSC (BSC classes grouped as 1+2, 3, 4, 5+6+7), dietary PC: principal component derived from clr-transformed percent energy dietary data, ASD PGS: polygenic score for ASD, ASD-ADHD-TS PGS: polygenic score from the ASD-ADHD-TS GWAS summary statistics, IQ-DQ composite score aggregated from MSEL non-verbal score, WISC-IV composite score, NIH Toolbox (age-corrected IQ); CSHQ: Children’s Sleep Habits Questionnaire raw score; ADOS2/G comparison score: composite of Autism Diagnostic Observation Schedule version 2 (ADOS-2) and -Generic version (ADOS-G) comparison score; ADOS2/G RRB and social affect CSS scores: ADOS-2 and ADOS-G Calibrated Severity Scores for RRB and social affect domains; SRS T-score: Social Responsiveness T-score; SSP sensor raw score: score from the “sensory” domain of the Short Sensory Profile. an = 96 UNR children with dietary PCs.
Results
Study characteristics
A total of 247 children (aged 2–17) participated in this study, including 99 diagnosed with ASD (“ASD”), 51 paired siblings without a diagnosis (“SIB”), and a total of 97 unrelated (including to each other) undiagnosed children without a diagnosis (“UNR”) (Figure 1A). Participants in the ASD and SIB groups came exclusively from the AAB, whereas the UNR group was comprised of 48 children from the AAB and 49 from QTAB. Two UNR participants had incomplete dietary data and were included when dietary variables were not analyzed. The AAB and QTAB stool samples were contemporaneously collected using an identical protocol and were receipted and processed by the same specialist human sample processing unit. The UNR_AAB participants (mean age 6.3, SD = 3.0) were on average younger than the UNR_QTAB participants (mean age 12.1, SD = 1.0), but after combining these sub-groups, the three study groups (ASD, SIB, UNR) were matched by age (ANOVA p = 0.2) (STAR Methods; Document S1). ASD studies commonly have a male sex bias, and after attempting to match groups by age and sex, the final percentage of males per group was ASD 73%, SIB 57%, UNR 54%. The AAB participants who provided stool samples for this study were generally representative of the wider AAB (Table 1). The ASD group on average had a lower intelligence quotient (IQ)-developmental quotient (DQ) composite score (IQ-DQ) (ASD mean = 85, SD = 24; SIB mean = 100, SD = 17; UNR mean = 100, SD = 15). As a measure of degree of ASD traits, the mean Autism Diagnostic Observation Schedule-2/G (ADOS2/G) comparison score in the ASD group was 6.7 (SD = 2.1). Full demographic details and comparisons are provided in Table 1, further detail on participant recruitment and the QTAB participants are provided in STAR Methods, and an overview of the phenotypic data used in this study is provided in Figure 1C.
Negligible variance in ASD diagnostic status is associated with the microbiome compared to age, stool, and dietary traits
First, we estimated the proportion of phenotypic variance associated with the microbiome for ASD diagnosis and a variety of other traits (Figure 1C), including neurodevelopmental traits (IQ-DQ, Children’s Sleep Habits Questionnaire [CHSQ] raw score), phenotypes with intuitive relationships with the microbiome (dietary PCs [Figure 1B], dietary diversity, stool consistency; see STAR Methods), and CD4+ T cell proportion, adjusting for covariates (described in Figure 2 legend). This variance estimate—the microbiome-association-index (or “b2”) (Rothschild et al., 2020)—is analogous to heritability (h2) from genetic analyses and provides an upper limit for predictive ability under the assumption of additivity. Whereas h2 reflects (direct or indirect) causality, b2 may reflect cause or consequence of trait variation. This analysis involves calculating an “omics relatedness matrix” (ORM) from microbiome features (taxa with a focus on the species-level and various hierarchies of microbial genes) between each pair of individuals and regressing against the trait in a linear mixed model framework (Zhang et al., 2019).
01231-9&id=gr2.jpg){kind=link}
We created ORMs from several measures of microbiome composition (Table S1) at the level of species and microbial genes (open-reading frames from the Microba Genes [MGENES] database, referred to in this section as “microbial genes”; Enzyme Commission classification [focusing on level 4 of this hierarchy] [Bairoch, 2000]; Transporter Classification Database [TCDB] [Saier et al., 2016]; MetaCyc pathways [Caspi et al., 2020]) and stratifying into common and rare features (Table S1; threshold based on median count > 0 for each feature). We examined the data distribution (Figure S1) and checked that the centered-log-ratio (clr) transformation was appropriate (Figure S2).
01231-9&id=figs1.jpg){kind=link}
01231-9&id=figs2.jpg){kind=link}
01231-9&id=figs3.jpg){kind=link}
We first analyzed age and BMI as benchmarking traits. Our results (in children) were consistent with previous species-level b2 estimates of 28% and 11% in adults, respectively (Rothschild et al., 2020): common species (n = 96) provided b2 estimates of 33% for age (SE = 8%, p∼0, false discovery rate [FDR]-significant; covariates: sex) (Figure 2; Table S1), and b2 = 12% for BMI (covariates: age, sex; SE = 7%, p = 3.5e-2, not FDR-significant). Notably, gene-level ORMs were associated with greater variance for both age (b2>99%, SE = 13%–17%, p ∼ 0, based on either n = 251,617 common genes with median count >0 or n = 1,742,727 rare genes with median count = 0) (Figure 2; Table S1) and BMI (rare genes b2 = 46%, SE = 22%, p = 8.4e-3, FDR-significant), suggesting that taxonomic and functional microbiome measures capture different elements of phenotypic variance.In notable contrast to the results for age, species- and gene-level b2 estimates for ASD diagnosis were weak and non-significant (maximum: rare genes b2 = 7%, SE = 16%, p = 0.33; covariates: age, sex, dietary PC1-3), as were those for IQ-DQ (common species b2 = 7%, SE = 13%, p = 0.39; covariates: age, sex), sleep problems measured by CSHQ raw score (common species b2 = 10%, SE = 9%, p = 0.17; covariates: age, sex) (Figure 2; Table S1), and clr-transformed CD4+ T cell proportion (b2∼ 0, SE = 0.06, p = 0.50; covariates: age, sex, dietary PC1-3) (Quantification and statistical analysis).Unlike these neurodevelopmental and immune traits, we observed strong FDR-significant b2 estimates for both stool consistency (covariates: age, sex, group, dietary PC1-3; rare species b2 = 41%, SE = 11%, p = 8.7e-6; rare genes b2 = 64%, SE = 20%, p = 2.5e-5) and dietary PC1 (rare genes: b2 = 48%, SE = 15%, p = 3.8e-4; covariates: age, sex, participant group) (Figure 2; Table S1; Quantification and statistical analysis).To explore the relative impact of related microbiome measures on b2 estimation as a sensitivity analysis, we explored fitting multiple ORMs. Fitting combinations of ORMs built from (1) taxonomic and functional datasets, (2) common and rare subsets of features, and (3) fitting multiple hierarchies (e.g., species, genus, and family data) all increased b2 estimates (Quantification and statistical analysis; Figure S3; Table S1).Overall, our results were robust to multiple sensitivity analyses: without covariates (Figure 2, open circles), a smaller age range, using the MetaPhlAn2 (Truong et al., 2015) taxonomic profiling pipeline, excluding those with current antibiotics usage, and excluding siblings (Methods S1). We also investigated non-additive models that estimate the predictive ability of the microbiome, again finding that the microbiome is predictive of age, but not ASD diagnosis (Quantification and statistical analysis).
Differentially abundant taxa and genes implicate Romboutsia timonensis
We next looked for microbial markers of ASD diagnosis, testing for differential abundance of 607 species, 297 genera, 38 orders, and 15 phyla. We used Analysis of Composition of Microbiomes (ANCOM)v2.1 (Mandal et al., 2015), a robust, non-parametric method that accounts for multiple testing and adequately controls the false positive rate (Weiss et al., 2017) (STAR Methods). Comparing ASD and the combined SIB and UNR groups (covariates: age, sex, and dietary PC1-3), only the species Romboutsia timonensis was significantly differentially abundant (lower abundance in ASD) at the conventional detection threshold > 0.7 (Figure 3A; Table S2). The results were consistent across extensive sensitivity analyses (Figures 3B and S4A–S4I; Methods S1), also finding that reduced Erysipelatoclostridium sp003024675 was often the next-most differentially abundant taxa (Table S2). In differential-presence testing (Fisher’s exact test, ASD versus UNR), the same two taxa (R. timonensis p = 3.9e-4, 56/99 ASD versus 78/97 UNR; E. sp003024675 p = 1.5e-4, 6/99 ASD versus 25/97 UNR) were also the top-ranked species, though neither survived FDR correction (Table S2). In permutation testing (n = 1000 random shuffles of diagnostic labels for each sample), these taxa were significantly differentially abundant (p ≤ 0.001) when compared to the empirical distributions for both ANCOM and Fisher’s exact tests (Table S2), adding further evidence that these findings are robust.
01231-9&id=gr3.jpg){kind=link}
01231-9&id=figs4.jpg){kind=link}
Notably, we failed to replicate previously reported ASD-gut microbiome associations (Ho et al., 2020) with the genus Prevotella, phylum Firmicutes, Clostridiales clusters, and species of Bifidobacterium (Table S2). However, we note that R. timonensis (family Peptostreptococcaceae, order Clostridiales, class Clostridia, phylum Firmicutes A) and E. sp003024675 (family Erysipelatoclostridiaceae, order Erysipelotrichales, class Bacilli, phylum Firmicutes) are members of these phylogenetic groups. Poor replication may be related to (1) prior microbiome studies being underpowered and prone to sampling biases (Rothschild et al., 2020); (2) technical differences between metagenomics and 16S rRNA sequencing—the former providing more detailed taxonomic resolution, which is relevant as R. timonensis was only recently isolated in human gut in the setting of anemia following bariatric surgery (Ricaboni et al., 2016); and (3) different statistical methods between studies, with variable use of adequate adjustment for multiple-testing and/or confounders.Unlike 16S sequencing studies (which have dominated the ASD microbiome literature), metagenomics sequencing permits functional insights. Thus, we looked for microbial genes associated with ASD diagnosis. We focused on relative abundances of metagenome annotations to MetaCyc groups, MetaCyc pathways, and Enzyme Commission (EC) gene families in order of increasing resolution (as opposed to the entire gene set due to computational and multiple-testing burden), finding no significant associations (covariates: age, sex, dietary PC1-3, Figures S4J–S4L) (Table S2).Next, we sought to identify specific genes or pathways from R. timonensis underlying the ASD-associated signal. We tested for ASD versus SIB+UNR differential abundance for MetaCyc groups, MetaCyc pathways, EC gene families, and specific genes (n = 4,950 genes with >10 non-zero values across samples) directly mapping to R. timonensis in our dataset (covariates: age, sex, dietary PC1-3) (Figures S4M–S4O). We identified six differentially abundant genes with a detection threshold > 0.7 (Figures 3C and 3D; Table S2), one of which overlapped with the EC gene set, whereas there were no associations with MetaCyc groups or pathways (Figures S4M–S4O; Table S2). Consistent with the species-level direction of effect, all significantly differentially abundant genes had reduced clr-transformed abundance in the ASD group (Figure 3D). Their functions included metabolism of amino acids (L-glutamine, L-lysine, L-methionine, and L-threonine), purines and pyrimidines, carbohydrates (galactose), as well as bacterial spore germination and dsDNA digestion (Quantification and statistical analysis; Table S2). We note that these results represent potential transcription of microbial genes and that metatranscriptomics data would be needed to evaluate actual expression.We performed species-level differential abundance analysis for IQ-DQ composite score, finding that lower Bifidobacterium sp002742445 passed detection threshold > 0.7 in the ASD versus SIB+UNR comparison but only passed detection threshold > 0.6 in the ASD versus UNR sensitivity analysis (Figures 3E and 3F; Table S2; Quantification and statistical analysis). There were no virome associations in the ASD versus SIB+UNR analysis (n = 200 taxa) (Table S2).
Dietary diversity mediates ASD-microbiome associations
Whereas ASD-associated signals were scarce in the metagenomics data, there were consistent associations with diet: in the variance component analysis, the dietary ORM was strongly associated with ASD (R2 = 14%, SE = 7%, p = 2.2e-5, FDR-significant) (Figure 2), and we observed significantly lower dietary PC3 in ASD (suggesting reduced meat intake as shown in Figure 1B) compared to SIB and UNR groups after adjusting for age and sex (Figure S5). Given that children on the spectrum favor less-diverse diets (Panossian et al., 2020; Schreck and Williams, 2006), we explored the effect of dietary diversity on the microbiome and gastrointestinal symptoms (via stool consistency).
01231-9&id=figs5.jpg){kind=link}
The ASD group had significantly less-diverse diet—estimated using Shannon index to measure dietary alpha-diversity from 123 food-level variables (see STAR Methods)—than both SIB and UNR groups (one-way ANOVA p = 1.3e-7, FDR-significant), including when adjusting for age and sex (one-way ANOVA p = 1.5e-6, FDR-significant) and energy intake (ANOVA p = 9.0e-7, FDR-significant) (Figure 4A; Table S3). They also had lower dietary quality, as measured using the validated Australian Recommended Food Score (ARFS) for children and adolescents (ANOVA p = 7.9e-4, FDR-significant) (Figure S6; Table S3; Methods S1) (Burrows et al., 2014; Marshall et al., 2012). These robust dietary results starkly contrasted with negligible direct association between ASD diagnosis and taxonomic alpha-diversity irrespective of measure and covariate use (age, sex, and dietary PC1-3) or beta-diversity (weighted Unifrac permutational multivariate analysis of variance [PERMANOVA] p = 0.20 with same covariates, PERMDISP2 p = 0.85) (Figure 5; Quantification and statistical analysis).
01231-9&id=gr4.jpg){kind=link}
01231-9&id=figs6.jpg){kind=link}
01231-9&id=gr5.jpg){kind=link}
Pursuing this link between ASD and dietary diversity, we hypothesized that taxonomic diversity may be a downstream consequence of diet rather than being directly associated with ASD diagnosis. Consistent with this hypothesis, there was a significant positive correlation between dietary and taxonomic diversity (Pearson r = 0.25, p = 6.3e-5, FDR-significant) (Figure 4C; Table S3), and in reciprocal regression analyses, dietary and taxonomic diversity (covariates: age, sex, stool consistency, and group) were significant predictors of each other (b = 0.02, p = 3.7e-2 and b = 1.1, p = 3.7e-2, respectively) (Figures 4D, 4E, S7A, and S7B with additional covariates). Furthermore, the largest effects in the dietary diversity regression were from group (UNR: b = 0.035, p = 3.0e-6; SIB: b = 0.021, p = 1.4e-2), whereas taxonomic diversity was not associated with group (Figures 4B and 4E). This suggests that ASD-associated dietary restrictedness (but not diagnosis itself) is associated with reduced microbiome diversity. These effects were robust to sensitivity analyses in which energy intake (kJ) was fitted as an additional covariate (Figure S7) and when substituting dietary diversity for dietary quality measured using the ARFS (Figure S6; Table S3; Methods S1).
01231-9&id=figs7.jpg){kind=link}
Next, we explored relationships between dietary and taxonomic diversity and stool consistency. We replicated a previously reported (Hadizadeh et al., 2017; Vandeputte et al., 2016; Zhernakova et al., 2016) inverse relationship between stool consistency and taxonomic diversity (b = −0.41 p = 3.7e-3 without covariates, FDR-significant) (Table S3) that was robust to covariates (Figure S8C; Table S3). We also identified a nominally significant association of stool consistency with dietary diversity (b = −2.34, p = 3.7e-2 in the model without covariates), although this did not survive covariate adjustment (Figure S8B; Table S3). Notably, the taxonomic diversity model (without covariates) explained greater variance (R2 = 3.2%) in stool consistency than the dietary diversity analysis (R2 = 1.5%), and in models of stool consistency with both taxonomic and dietary diversity fitted as explanatory variables, only taxonomic diversity was significant (b = −0.36, p = 1.5e-2) (Figure S8A) when covariates were included (Figure 4F). Overall, this suggests that looser stool consistency (higher rBSC score) is proximally related to reduced taxonomic diversity, which is downstream of reduced dietary diversity. This mechanism may explain the reported relationship between increased gastrointestinal issues and increased repetitive behaviors (Chakraborty et al., 2020).
01231-9&id=figs8.jpg){kind=link}
Behavior and preferences are upstream of reduced dietary and taxonomic diversity
We investigated whether behavioral factors diagnostic of ASD are upstream of restricted diet and reduced dietary and taxonomic diversity. To achieve this, we leveraged both psychometric measures and polygenic scores from human genotyping data (for ASD and other phenotypes), the latter using genome-wide SNP data generated from blood samples in the AAB (Yap et al., 2021) and QTAB participants (see STAR Methods). Whereas previous sections discussed the contributions of microbial genes, this section explores the relationship between human polygenic scores and dietary and taxonomic diversity.First, we confirmed the ASD-dietary diversity association through analysis of continuous autism-spectrum measures. We identified an inverse association between ASD polygenic score in this cohort (Yap et al., 2021) and dietary diversity (b = −1.0e2, p = 1.2e-2, FDR-significant), but not taxonomic diversity (b = −4.4e-2, p = 0.17) (Figure 4G). We note that the ASD polygenic score itself was insufficiently powered to predict ASD diagnosis in this dataset (one-way ANOVA p = 0.41; for results in the entire AAB refer to Yap et al., 2021). Phenotypically, we observed negative associations between dietary diversity and two quantitative measures of degree of ASD features: ADOS2/G comparison score (b = −8.6e-3, p = 3.1e-3, FDR-significant; n = 99 ASD) (Figures 6A and 6B) and Social Responsiveness Scale t-score (b = −6.4e-4, p = 7.8e-2, marginally significant; n = 97 AAB children: 10 ASD and 87 SIB/UNR) (Figures 6C and 6D).
01231-9&id=gr6.jpg){kind=link}
Second, we hypothesized that repetitive and restrictive behaviors and interests (one of two diagnostic domains for ASD) may underlie a restricted diet, upstream of microbiome changes. Phenotypically, we observed FDR-significant negative association between higher combined (Lord et al., 2012) ADOS-2/G restricted and repetitive behavior (RRB) Calibrated Severity Scores and dietary diversity (without covariates: b = −7.8e-3, p = 3.8e-3; with covariates age and sex: b = −6.4e-3, p = 1.8e-2; both FDR-significant) (Table S3) and nominally-significant negative association with taxonomic diversity (without covariates: b = −4.3e-2, p = 2.4e-2; with covariates age and sex: b = −3.4e-2, p = 7.6e-2) (n = 99 ASD group only) (Figure 4H; Table S3). We then leveraged genome-wide association study (GWAS) summary statistics from a cross-trait analysis (Yang et al., 2021) of ASD, attention-deficit hyperactivity disorder, and Tourette Syndrome (hereafter called ASD-ADHD-TS) to generate polygenic scores (PGS) for restrictive-repetitive behaviors. We confirmed that ASD-ADHD-TS PGS correlated with ADOS-2/G RRB scores in the full AAB (r = 0.09, p = 0.01, n = 867) and so represents a genetic proxy. We found marginal association between ASD-ADHD-TS PGS and reduced dietary diversity (b = −7.2e-3, p = 0.10) but not with taxonomic diversity (b = −4.7e-2, p = 0.18) (Figures 6E and 6F). We compared these results to ADOS2/G social affect Calibrated Severity Scores, the other major domain of ASD diagnosis, noting that ADOS2/G RRB and social affect scores were significantly correlated (r = 0.29, p = 3.7e-3, n = 99 participants). We again found FDR-significant associations between ADOS2/G social affect and dietary diversity (without covariates: b = −8.2e-3, p = 1.0e-2; with covariates age and sex: b = −7.5e-3, p = 1.5e-2 both FDR-significant) (Figure 6G; Table S3) and not with taxonomic diversity (without covariates: b = −3.3e-2, p = 0.15; with covariates age and sex: b = −2.9e-2, p = 0.19) (Figure 6H; Table S3). Importantly, these associations were weaker than with the ADOS2/G RRB scores. Thus, ASD-associated restricted and repetitive behaviors appear to have stronger relationships with dietary diversity than social affect.Third, on the basis that sensory sensitivity may also underlie restricted dietary preferences (Cermak et al., 2010), we explored relationships with Short Sensory Profile raw sensory score in a small ASD-only subset of the data for which this instrument was completed. We found marginal associations with both dietary (b = −9.5e-4, p = 6.9e-2) and taxonomic diversity (b = −6.8e-2, p = 8.6e-2) (n = 91; Figure 4I). Notably, the Short Sensory Profile raw sensory score was not correlated with ADOS2/G RRB Calibrated Severity Scores in this dataset (r = 0.05, p = 0.64).In contrast, we found no evidence for hypothesized links between dietary and taxonomic diversity measures and ASD-associated traits such as neuroticism as a proxy for anxiety (Kim et al., 2018; Yang et al., 2019) (Figures 6K and 6L; Table S3; Quantification and statistical analysis).Overall, these data suggest that ASD-associated preferences and behaviors lead to reduced dietary diversity, which mediates weak ASD-microbiome relationships (Figure 4J). Notably, all psychometric measures had more significant correlations with dietary diversity than taxonomic diversity (Figures 4G–4I and 6A–6H). However, we cannot rule out the possibility that these downstream microbiome effects could also feed-back and influence behavior.
Discussion
In this large ASD stool metagenomics study in which confounders were carefully considered, we found negligible evidence for direct associations between the stool microbiome and ASD diagnostic status, which was also the case for other neurodevelopmental traits (e.g., IQ-DQ, sleep problems). For ASD, there was limited evidence for associations with taxonomic diversity or microbiome-association index (b2; Figure 2), and only one differentially abundant species was robustly identified (Figure 3). These results were striking when compared to strong associations of microbiome composition with age, diet, and stool consistency (Figure 2). Importantly, we failed to replicate previously reported ASD-microbiome associations from human studies. Instead, we found evidence linking behaviors associated with the autism spectrum (e.g., repetitive-restricted behaviors or interests, sensory preferences, and social affect) to reduced dietary diversity, which, in turn, was associated with reduced microbiome diversity and looser stool consistency (Figure 4J). This putative model challenges suggestions from animal studies that the microbiome may be causally related to ASD-related traits (Buffington et al., 2016; Hsiao et al., 2013; Sharon et al., 2019). Our findings also stand at odds to the proliferation of experimental interventions and early clinical trials that propose to “treat” ASD by targeting the microbiome (Kang et al., 2019; Kang et al., 2017).In contrast to measures of microbiome composition, ASD was robustly and significantly linked to dietary variables, irrespective of covariates (Table S3). We found (1) that significant variance in ASD diagnosis was associated with diet but not the microbiome in the b2 analysis (Figure 2), (2) reduced meat intake in the ASD group (Figure S5), and (3) reduced dietary diversity in the ASD group despite significantly higher variance in dietary diversity (Figure 4A), which is consistent with the dietetics literature (Panossian et al., 2020; Schreck and Williams, 2006) and some smaller ASD microbiome studies with dietary data (Berding and Donovan, 2018).One rationale for the interest in the ASD microbiome is the frequent co-occurrence of gastrointestinal complaints (Chaidez et al., 2014; Kohane et al., 2012; McElhanon et al., 2014). In the absence of complete gastrointestinal complaint reporting, we analyzed stool consistency scores, with the caveat that it is unclear how this single-time point data reflects chronic conditions. Stool consistency appeared to be more proximal to taxonomic than dietary diversity, although we acknowledge that it is difficult to distinguish between a top-down (i.e., dietary and taxonomic diversity influencing downstream stool consistency) versus bottom-up (i.e., stool consistency being an upstream proxy) relationship. For the former, dietary restrictedness could plausibly affect gut ecology and taxonomic diversity, which in turn affects stool consistency. In relation to a bottom-up model, looser stool may indicate underlying food allergies or intolerances, which may be associated with deliberate (parental) dietary restriction to identify causative agents. Additionally, looser stool consistency reflects reduced gastrointestinal transit time and reduced colonic water reuptake (Vandeputte et al., 2016), which affects taxonomic diversity. As the narrow-sense heritability of gastrointestinal conditions that affect stool consistency (e.g., irritable bowel syndrome) are small (Wu et al., 2021), environmental contributions likely predominate over genetics (Rothschild et al., 2018).Our results have important implications for understanding the role of the gut microbiome in ASD and other psychiatric traits. First, in relation to medical care, food selectivity among children on the autism spectrum is an important clinical concern. Food selectivity is related to avoidant/restrictive food intake disorder (ARFID; which is likely underdiagnosed despite affecting over 20% of autistic children [Koomar et al., 2021]) and can cause nutritional deficiencies among autistic children (Zimmer et al., 2012) to the extent that hospitalization and invasive measures such as enteral feeding are required (Tang et al., 2011). Our results also suggest that dietary quality is poorer in children on the spectrum (Methods S1). Given that elevated microbial diversity is robustly associated with improved health outcomes (Valdes et al., 2018), the association of ASD with poorer dietary quality and reduced dietary and taxonomic diversity underlines the importance of dietary and nutritional interventions in this population. Second, our results have implications for the interpretation of cause and effect in relation to diet in microbiome analyses in psychiatric conditions. There is growing interest in the contribution of diet and the microbiome to psychiatric traits (e.g., depression [Dash et al., 2015; Molendijk et al., 2018]), but our results emphasize the need to consider the (arguably more intuitive) impact of behavior on the microbiome (Jacka et al., 2015). These results add to other reports of diet driving microbiome associations with health (Claesson et al., 2012).For future microbiome studies, we emphasize the importance of collecting detailed dietary data (recent examples [Asnicar et al., 2021; Wang et al., 2021]), particularly for ASD and other neuropsychiatric traits with plausible co-variation of diet with diagnosis or treatment. We advocate for larger sample sizes to ensure that results are robust to sampling effects and to identify subtler microbiome associations. We also recommend higher-resolution metagenomics technology and expanded databases since more granular taxonomic measures of microbiome composition were more sensitive (Table S1), gene-level ORMs explained more variance for some traits (Table S1), power to detect associations was weaker with the MetaPhlAn2/NCBI pipeline (Methods S1), and because taxonomic and functional datasets may capture complementary aspects of the microbiome (Figures S1 and S3).In conclusion, we found negligible direct associations between ASD and the gut microbiome in contrast to strong associations with other phenotypes such as age, dietary variables, and stool consistency. Instead, we find evidence that restricted dietary diversity and poorer quality—which is associated with specific ASD features such as restrictive-repetitive behaviors—is a significant mediator of taxonomic diversity, and in turn, stool consistency. Our results are consistent with an upstream role for ASD-related behaviors and dietary preferences on the gut microbiome and are contrary to claims of the microbiome having a major (or causal) role in ASD.
Limitations of the study
First, this study did not have a longitudinal design, so we cannot rule out microbiome contributions prior to ASD diagnosis. Second, although this is to our knowledge the largest metagenomics study of the ASD stool microbiome to date, there may nonetheless be sampling biases that require larger studies to overcome (Rothschild et al., 2020). Third, this study used stool samples as a gut microbiome proxy, which may not accurately represent the more difficult-to-access mucosal microbiome (Shanahan et al., 2016). Fourth, data on antibiotic intake in this cohort were not systematically collected and so could not be rigorously accounted for other than through exclusion in sensitivity analyses. Fifth, the gold-standard differential abundance analysis relied on per-feature tests that do not reflect the interactions and non-independence that occurs within an ecological or metabolic context. Finally, we await the emergence of datasets with comparable study design, consideration of confounders, and depth of phenotypic and metagenomics data for replication of these results.
STAR★Methods
Key resources table
Resource availability
Lead contact
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Jacob Gratten ( jacob.gratten@mater.uq.edu.au ).
Materials availability
This study did not generate new unique reagents.
Experimental model and subject details
This was a human study. The 247 participants in this study came from two datasets: the Australian Autism Biobank (AAB) (Alvares et al., 2018) and the Queensland Twin Adolescent Brain (QTAB) project. The AAB participants include children with an ASD diagnosis (recruited from autism clinics and research centers across Australia’s four largest cities: Sydney, Melbourne, Brisbane, Perth; no exclusion criteria; mean age = 8.7, SD = 3.8, 73% male), and their siblings (“SIB”) without a diagnosis (mean age = 8.0, SD = 4.3, 57% male) (Alvares et al., 2018). The group of unrelated children (“UNR”) without a diagnosis were recruited by the AAB (recruited from the community; exclusion criteria: having an ASD diagnosis) and QTAB (typically developing children recruited from the community) (AAB+QTAB UNR group mean age = 9.2, SD = 3.7, 54% male). Stool samples from AAB and QTAB participants were collected using the same protocol, in a contemporaneous time frame, and were receipted and processed by the same University of Queensland Human Studies Unit, prior to being randomized for DNA isolation, library preparation and DNA sequencing by Microba Life Sciences. Overall, the sample included n = 51 sibling-pairs overlapping with a case-control design comprising n = 99 cases and n = 97 unrelated controls (Table 1,Figure 1a). A description and comparison of demographic data from participant groups is provided in Table 1. The AAB subset included in this metagenomics study was representative of the wider AAB (Table 1).Sample selection for this study was constrained to AAB and QTAB participants who provided a stool sample and completed a dietary questionnaire. Within these constraints, we attempted to match by sex and age. Specifically, we selected UNR_QTAB participants to ameliorate a disparity in sex (chi-square statistic = 8.3, p = 1.5e-2) and age (ANOVA p = 2.0e-3) between the UNR_AAB subset (n = 48 participants), who were younger (mean age 6.3, SD = 3.0) with equal sex balance (50%), and the ASD (AAB) group (n = 99), who were older (mean age = 8.7, SD = 3.8) with a male bias (73%). In the final sample, there was a persisting (albeit reduced) male bias (chi-square statistic = 7.7, p = 2.2e-2), consistent with the elevated male:female ratio in ASD, but no significant difference in age between the groups (ANOVA p = 0.2) (Table 1). For this reason, we statistically adjusted for sex in all analyses. Overall, including the QTAB participants in the study increased power, and ameliorated age (and to a lesser degree, sex) differences that are an important modulator of the microbiome, with minimal likelihood of collection bias.
Ethics approval and consent to participate
All families provided informed consented to be included within this study.• NSW: Sydney Children’s Hospital Network HREC, approval number HREC/14/SCHN/269.• QLD: Mater Health Services HREC, approval number HREC/14/MHS/212; the University of Queensland, approval number 2014001079; QTAB Project: Children’s Health Queensland HREC, approval number HREC/16/QRCH/270; The University of Queensland, approval number 2016001784/ HREC/16/QRCH/270• VIC: La Trobe University, approval number HEC16/104• WA: Princess Margaret Hospital for Children approval number 2014029EP; the University of Western Australia approval number RA/4/1/8184
Method details
Phenotype data
Dietary data
Dietary data were collected in both AAB and QTAB cohorts – predominantly reported by parents – using the Australian Eating Survey (AES; toddler and children’s versions) (Collins et al., 2013; Watson et al., 2009), which has been validated in the Australian population. Food-level intake data were available for n = 245 of the 247 participants, and percent energy (pe) data were available for n = 246. The AES records frequencies of intake for 123 different foods, from which derived variables are generated, including pe from each of 13 core (vegetables, fruit, meat, alternative proteins, grains, dairy) and non-core (sweet drinks, packed snacks, confectionery, baked products, takeaway, condiments, fatty meats) food groups; macronutrients (various carbohydrates, fats and proteins); micronutrients (various vitamins and minerals), and the Australian Recommended Food Score (ARFS).We primarily used the dietary data in two ways (Figure 1c).1. We used the food-level input to measure dietary diversity (n = 245) using Shannon index – the same measure of alpha-diversity used in the microbiome analyses. The AES records the frequency of intake for each food item using ordinal variables, with the exact variable depending on the food item (e.g., for fruit intake: “Never,” “Once per week,” “2-4 per week,” “5-6 per week,” “1 per day,” “2 or more per day”; whereas for hamburgers: “Never,” “Less than 1 per month,” “1-3 per month,” “1 per week,” “2-4 per week”). We encoded these N categories on an integer ordinal scale, where the least-frequently consumed category was encoded as 1, and the most-frequently consumed category was encoded as N. We then calculated dietary diversity from this integer dataset using the Shannon index. For the variance component analysis, we input the integer dataset to generate an ORM (correlation matrix between all pairs of individuals) for the OREML analysis. To further test the validity of our dietary diversity construct, we also tested the Australian Recommended Food Score for children and adolescents (Burrows et al., 2014; Marshall et al., 2012), which is a validated measure of dietary quality (Methods S1).2. We calculated principal components from the percentage energy data (n = 246 participants; referred to as dietary PCs; STAR Methods, Methods details), to capture salient dietary features that may affect the microbiome, given that a strong relationship has been identified by others (Bokulich et al., 2016; Rothschild et al., 2018).On the basis of their loadings onto the dietary items (Figure 1b), the first 3 dietary PCs may be interpreted as follows:• PC1: a high value is associated with a diet high in plant-based foods (vegetables, fruit, alternative proteins) and low in non-meat non-core foods (sweet drinks, packed snacks, confectionery, baked products, takeaway, fatty meats)• PC2: a high value is associated with a diet high in dairy products and low in grains and takeaway• PC3: a high value is associated with a diet high in meat (including fatty meats) and low in grains and dairy.An RMarkdown document showing exploratory data analysis of the dietary data is provided at https://zenodo.org/record/5558047.
Bristol Stool Chart
The Bristol Stool Chart (BSC) was completed by each AAB and QTAB participant in the study. The BSC provides a pictorial scale of stool consistency numbered 1 to 7, whereby lower values indicate dryer stool (e.g., constipation) and higher numbers denote watery stool (e.g., diarrhea). As we had limited representation at the extremes of this distribution, we regrouped scores of 1+2 and 5+6+7 together (Figure 1c). This resulted in a 4-point scale with > 50 samples per group. We treated this regrouped BSC variable (rBSC) as a continuous variable.
Neurodevelopmental phenotypes
Our analyses included multiple neurodevelopmental phenotypes. ASD-related psychometric measures included either ADOS-2 or ADOS-G repetitive and restricted interests sub-scores and social affect sub-scores (combined as per (Lord et al., 2012), and with domains converted to the Calibrated Severity Score to allow comparison across modules (Hus et al., 2014; Hus and Lord, 2014)), Social Responsiveness Scale (Constantino, 2002) t-score (undiagnosed AAB children), and the Short Sensory Profile (McIntosh et al., 1999) raw Sensory score. Among non-ASD-related traits, we generated a composite IQ (from WISC-IV (Wechsler, 2003) for older children in the AAB, and the NIH Toolbox age-corrected full-scale score (Akshoomoff et al., 2013) for the QTAB study) and non-verbal developmental quotient score (from MSEL (Mullen, 1995) for younger children in the AAB). These questionnaires were aggregated to capture a larger proportion of the dataset, and to provide a proxy for intellectual and developmental delay, which we refer to hereafter as IQ-DQ. We combined these scores to improve power, and all these measures were similar in that they were normed to approximately conform to the expected mean = 100 and SD = 15 in the population. We also looked for association with sleep disturbances in the AAB subset, using the Children’s Sleep Habits Questionnaire (CSHQ) (Owens et al., 2000).
Polygenic scores (PGS) from human genotyping
We calculated PGS for ASD (Grove et al., 2019), ASD-ADHD-TS cross-disorder analysis (Yang et al., 2021), and neuroticism (Nagel et al., 2018) in the AAB and QTAB datasets, using the pipeline described in (Yap et al., 2021). Briefly, we filtered for genome-wide association study SNPs that were in the HapMap3 reference, and were common across both datasets, then input all of these filtered SNP summary statistics into SBayesR (Lloyd-Jones et al., 2019) to re-weight the SNP effect sizes (settings: –pi 0.95, 0.02, 0.02, 0.01; –gamma 0, 0.01, 0.1, 1; –chain-length 10000; –burn-in 2000; –out-freq 10, and using the –exclude-mhc flag). We then multiplied the best guess genotypes in the target sample (i.e., AAB individuals and UKB controls) by the re-weighted effect sizes, using the PLINK (Chang et al., 2015; Purcell and Chang, 2015) –score function.
Cell-type proportions
We calculated cell-type proportions for the subset of participants for whom we had also generated DNA methylation data. For this, we applied the Houseman algorithm implemented in meffil (Min et al., 2018), using the default blood reference (Reinius et al., 2012), which provides cell-type proportions for neutrophils, monocytes, B cells, CD4+ T cells, CD8+ T cells, NK cells, and eosinophils.
Sample collection and preparation
The same protocol was applied across both AAB and QTAB datasets.
Sample collection
As described in (Alvares et al., 2018), teaspoon-sized stool samples were collected by parents at home, either scraped from diapers or from a liner suspended in a toilet bowl, and suspended in 4mL RNAlater. Samples were brought to the clinic and shipped to the Institute for Molecular Bioscience at the University of Queensland the same day. In the vast majority of cases, stool samples were received within 2-3 days of collection (1-2 days from shipping). Time from shipping to processing was 12-72 h. Processing involved vigorous homogenization, before aliquoting and long-term storage at −80°C. All samples underwent only one freeze-thaw cycle, which was at the time of sequencing.
DNA extraction
DNA extraction was performed using a modified protocol, optimizing the initial mechanical lysis step, before following the manufacturer’s instructions for the QIAamp 96 PowerFecal QIAcube HT Kit (QIAGEN) on the QIACube HT automated extraction system. Prior to processing, samples were transferred to a 96-well plate and two washes with ice cold PBS performed as per the manufacturer’s instructions.
Library preparation
Libraries are prepared according to the manufacturer’s protocol using Nextera XT Library Preparation Kit (Illumina #FC-131-1096) with reduced reaction volume to allow processing in 384-well format. The libraries were indexed with NexteraXT v2 384 Index A-D (Illumina FC-131-2001-4). Resulting libraries were quantified and assessed for appropriate QC including fluorometric quantification and gel analysis.
Library pooling, QC, loading, and sequencing
Nextera XT libraries were pooled at equimolar amounts to create a sequencing pool. The pool was quantified and pool QC performed with gel analysis and fluorometric quantification. The library was prepared for sequencing on the Illumina NovaSeq6000 according to manufacturer’s instructions and sequenced with 2 × 150bp paired-end chemistry in the Microba laboratory. Pools were sequenced to a target depth of 3Gb per sample with a minimum of 2GB (approximately 7M – 16M paired-end reads).
Quantification and statistical analysis
Metagenomics datasets and QC
Quality control
Metagenomic sequencing data QC and processing was performed by Microba Life Sciences Limited. Paired-end DNA sequencing data were demultiplexed and adaptor trimmed using Illumina BaseSpace Bcl2fastq2 (v2.20), accepting one mismatch in index sequences. Reads were then quality trimmed and residual adaptors removed using the software Trimmomatic (v0.39) with the following parameters: -phred33 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 CROP:100000 HEADCROP:0 MINLEN:100. Human DNA was identified and removed by aligning reads to the human genome reference assembly 38 (GRCh38.p12, GCF_000001405) using Burrows-Wheeler Aligner (BWA) v0.7.17 (Li and Durbin, 2009) with default parameters except minimum seed set to -k 31. Alignments were further filtered using SAMtools v1.7 (Li and Durbin, 2009), with flags -ubh -f1 -F2304. Any remaining pairs where at least one read mapped to the human genome with > 95% identity over > 90% of the read length were flagged as human DNA and removed. Table S4.1 and Table S4.2 provide summary statistics of the QC process. There were between 4.8-24.8 million reads passing QC per sample (Table S4.2). All samples were then down-sampled to a standard depth of seven million read pairs to mitigate the effect of increased sensitivity in our association analyses. For n = 19 samples with fewer than 7 million reads, we retained all QC-ed reads, rather than excluding these samples, noting the ongoing debate on the need for rarefaction (McMurdie and Holmes, 2014).
Quantification of microbial abundances
Species profiles were obtained with the Microba Community Profiler v2.0.2 using the Microba Genome Database (MGDB) v2.0.0 as a reference database (Parks et al., 2021). Briefly, reads are mapped to MGDB, and the relative cellular abundance of species clusters with sufficient evidence to be considered reliable was estimated and reported. In this dataset, a total of 1,757 species were identified (964 specific to MGDB), including seven that were newly mined from this dataset. Of these, there were 1,735 bacteria, 12 archaea and 10 eukaryota.Quantification of gene and pathway abundance in the metagenomic samples was performed using the Microba Gene and Pathway Profiler (MGPP) v0.1.0 against the Microba Genes (MGENES) database v2.0.0. MGPP is a two-step process. In step one, all ORFs from all genomes in MGDB were clustered against UniRef90 release 2019/09 using 90% identity over 80% of read length using the tool MMSeqs2 Release 10-6d92c (Steinegger and Söding, 2017). Gene clusters were then annotated with the UniRef90 identifiers and linked to the Enzyme Commission (EC) and Transporter Classification Database (TCDB) annotations via the UniProt ID Mapping service. EC annotations were used to determine the encoding of MetaCyc (Caspi et al., 2020) pathways in each genome using enrichM [https://github.com/geronimp/enrichM] and pathways that were complete or near complete (completeness > 80%), were classified as encoded and stored for further analysis. In step two, all DNA sequencing read pairs that align by one or more base (in nucleotide space) to a gene sequence from any gene within a MGENES protein cluster were summed and tabulated. Pathway abundances were calculated by averaging the gene counts of each pathway present within all genomes of all species reported by MCP. In this dataset, a total of 5,165,783 genes were identified.
Viral species
We additionally profiled the virome (rarefaction to a constant 4.8 million reads) using MiCoP (LaPierre et al., 2019), as this method is optimized to call viruses and eukaryotes. As a reference dataset, MiCoP draws upon the NCBI’s RefSeq Viral database (Brister et al., 2015). We identified a total of 200 viral species in the dataset.
Transformations and filters
Metagenomics data (taxonomic and functional) are a form of compositional data, which violate assumptions of independence as proportionality imposes negative correlations within the dataset. Hence, we applied centered-log-ratio (clr) transformation (Aitchison, 1982) in the variance component analyses (offset = 0.001), whereas for differential-abundance analysis, additive-log-ratio transformation was performed upon the count data within the ANCOMv2.1 package (Mandal et al., 2015) (offset = 1). For the variance component and differential abundance analyses, we removed ultra-low prevalence features with < 10 non-zero values. This left 607 bacterial species (221 specific to the Microba database), 40 viral species, and 1,742,729 genes. However, we retained all features for Fisher’s exact test and diversity analyses.
Covariate choice
We performed preliminary analyses to identify covariates for inclusion in our various analyses. This process, which is described in greater depth below, led us to include age, sex and the top 3 principal components calculated using the derived variables of (clr-transformed) percentage energy from various food groups (“dietary PCs”). Together, these covariates explained 13.5% of variance in microbial taxonomic alpha-diversity (Shannon index) within the sample.
Demographic variables
Age and sex were included as baseline covariates. Month was initially included (as a factor variable) to account for seasonal variation in diet, but was subsequently dropped as it made minimal contribution to both measures of dietary intake and taxonomic alpha-diversity (Shannon index) (an additional 3.8% of variance comparing the age + sex model which explained 12.2% of variance and the age + sex + month model which explained 16.0% of variance). A seasonality variable that transformed month of dietary survey to a cosine curve made no contribution to the variance in taxonomic alpha-diversity (Shannon index).
Dietary principal components
We included dietary data from the Australian Eating Survey as covariates in the metagenomics analyses, specifically focusing on the percentage energy measures. This was because this metric best accounted for the confounding effect of energy intake (which is strongly correlated with age), and also had the highest completeness rates.To extract salient features of the dietary data, we generated principal components (dietary PCs) on the 13 percentage energy variables. As this measure is a form of proportional data – and therefore requires specific compositionally-aware methods – we performed a centered-log-ratio (clr) transformation before calculating the PCs. The importance of using the clr-transformation was suggested based on three lines of evidence. First, we found that age, sex and month covariates explained sizeable differences in variance for dietary PC1-3, depending on whether a clr-transformation was applied. Second, the correlation between the non-transformed and clr-transformed dietary PCs was low (r = 0.32), suggesting that compositionally-aware analysis may provide different results from a compositionally-unaware analysis. Third, the clr-transformed data explained a greater proportion of variance in microbial taxonomic alpha-diversity than the non-transformed data (percentage energy items explained 3.7% of variance with clr-transform, versus 2.5% of variance without clr-transform).The first 3 dietary PCs accounted for 43% of dietary variance. Including the first 3 clr-transformed PCs as covariates explained greater variance in microbial taxonomic alpha-diversity than all 13 percentage energy items (2.8% versus 2.5%), justifying the former’s inclusion as covariates in the metagenomics analysis. We note that the ASD group had lower dietary PC3 (suggesting reduced meat intake) after adjusting for age and sex (Figure S5).
Variance component analysis
Omics-relationship matrices (ORM) and OREML
We performed variance component analysis using the software package OSCA (Zhang et al., 2019). Briefly, this method estimates an n x n omics-relationship matrix (ORM) between each pair of n individuals based on p probes or variables (in this analysis, clr-transformed counts of taxonomic or functional variables from the metagenomics dataset, after removing variables with < 10 non-zero values); we used the –orm-alg 2 setting. Then, OREML was used to estimate the proportion of variance of a given phenotype (the dependent variable in the model) associated with the ORM, fitted as a random effect in a restricted maximum likelihood (REML) framework. Covariate choice depended on the focal phenotype, but universally included sex and age (except when age was the dependent phenotype), and in some cases, participant group or dietary PCs.We generated separate ORMs for common features (those with median count > 0) versus rare features (with median count = 0, but with ≥ 10 non-zero counts in the full sample) (Table S1.1). This was intended to approximate a mixture model when fitting these two ORMs in a multiple ORM framework (–multi-orm), motivated by the observation that there are some “core” taxa (approximately normally distributed) as well as “accessory” taxa (that are less prevalent and which have a zero-inflated distribution). Covariates per analysis are listed in the Figure 2 legend.We performed extensive sensitivity analyses for the variance component analyses: without covariates, restricting to AAB participants younger than 10 years of age, using a different taxonomic and functional profiling pipeline (MetaPhlAn2 (Truong et al., 2015) and HUMAnN2 (Franzosa et al., 2018)), and excluding participants with current antibiotic use. These analyses are detailed in Methods S1.
clr-transformation versus 0/1 coding for ORMs
Given that there was large spread in the distribution of rare species, we compared the effect of applying a clr-transformation to taxa counts versus encoding presence/absence (0/1) for ASD diagnosis, stool consistency (rBSC) and dietary PCs. We found that results were highly similar, supporting our choice to use the clr-transformed data (Figure S2).
CD4+ T cell proportion and variance analyses
CD4+ T cell proportions in blood were predicted for n = 151 participants from blood-derived DNA methylation data, using the meffil package. This was of interest because others have reported interplay between the microbiota and CD4+ regulatory T cells (Zheng et al., 2020) and because immune conditions commonly co-occur with ASD (Atladóttir et al., 2010; Sabourin et al., 2019; Vargason et al., 2019). To reduce multiple-testing burden, we initially only tested the effect of common genes (as this dataset typically had large effects), finding that there were no relationships (b2∼0, SE = 0.06, p = 0.50).
Dietary traits and variance analyses
Microbiome composition was associated with variance in dietary PC1, with b2 = 13% (SE = 7%, p = 0.03) when microbiome composition is based on common species, b2 = 8% (SE = 5%, p = 0.4) based on rare species, b2 = 48% (SE = 15%, p = 8.1e-5, FDR-significant) based on common genes and b2 = 22% (SE = 18%, p = 7.4e-2) based on rare genes (Figure 2, Table S1.2). Estimates from the microbiome ORMs were lower for dietary PC2-3 compared to PC1 (Figure 2, Table S1.2). However, the food ORM was consistently associated with the highest proportions of variance for dietary PC1-3 (b2∼60%), and dietary diversity (b2 = 52%; Figure 2, Table S1.2). Further results for the association between dietary PC1 and combinations of multiple ORMs are provided in the “Variance in traits explained by combinations of ORMs” section below.
Effect of age on age-associated dietary PCs
We noted that dietary PC1 (representing high percentage energy from a plant-based diet) was strongly associated with age. Hence, we tested the effect of removing age as a covariate in the OREML analysis on b2 estimates, finding that the b2 estimates were increased across multiple ORM datasets as including these covariates reduced the phenotypic variance (the denominator in the b2 estimate). Specifically, removing age as a covariate increased variance estimates from 49% to 55% for the common genes dataset, 21% to 38% for the rare genes dataset and from 13% to 18% for the common species dataset.
Effect of multiple ORMs on variance estimates
Benchmarking traits. We used dietary PC1 as an internal benchmarking trait because (1) dietary PC1 reflects a plant-based diet, which is known to be related to the gut microbiome (David et al., 2014) (Figure 1b), (2) it explained the most variance in the dietary PC analysis, and (3) variance in dietary PC1 was strongly associated with species-level microbiome composition (Figure 2, Table S1.2-S1.3).Combining ORMs within hierarchies. We evaluated the effect on variance explained in dietary PC1 after collapsing taxa into higher taxonomic levels (species versus genera versus family) and functional levels (for the Enzyme Commission dataset, comparing level 4 and level 3), taking the “common” variables (defined as median count > 0 in this dataset). In both of these taxonomic (species: b2 = 13%, SE = 8%, p = 0.03; genus: b2 = 14%, SE = 7%, p = 0.01; family: b2 = 4% SE = 4%, p = 0.15) and functional (Enzyme Commission level 4: b2 = 48% SE = 13%, p = 3.8e-3; Enzyme Commission level 3: b2 = 7% SE = 7%, p = 0.50) datasets, the most-granular data hierarchies tended to explain the greatest proportion of variance (Figure S3, Table S1.3), so we hereafter focus on these. We found b2 = 24% (SE = 9%) for dietary PC1 when combining each of the species, genus and family ORMs into the model (Figure S3, Table S1.3). These results also suggest that ORMs within taxonomic and functional hierarchies are not orthogonal.Calculating ORMs based on all features. We tested the effect of providing as input an ORM based on all species’ clr-transformed data (i.e., not stratifying ORM calculation into common versus rare features). We found that these estimates were highly similar to the multiple ORM analysis that included both common and rare species ORMs (e.g., for age, the “species_all” ORM explained b2 = 59% SE = 9%, whereas the common+rare species ORM (2-ORM model) explained b2 = 60% SE = 9%; whereas the individual common species ORM explained b2 = 33% SE = 8%, and the individual rare species ORM explained b2 = 53% SE = 9% further detail in Figure S3, Table S1.3). For ASD diagnosis, the “species_all” ORM explained b2 = 0% (SE = 8%), whereas the common+rare species ORM explained b2 = 1% (SE = 9%). This suggests that the ORM stratification is a valid approach.Common and rare features. We initially stratified common and rare taxa to capture core versus accessory variables, and to provide more granular insights into their relative contributions. We tested whether the combination of these ORMs (approximating a mixture model) could improve the variance explained for a given trait (in this case, dietary PC1), finding that the estimates were essentially equivalent (including common+rare species gave b2 = 23%, SE = 12%; versus common species alone b2 = 13%, SE = 8% and rare species alone b2 = 8%, SE = 9%). Interestingly, combining common and rare species ORMs provided a similar variance estimate for dietary PC1 compared to the combination of ORMs calculated using common features from emergent taxonomic hierarchies (common + rare species b2 = 23%, SE = 12% versus species + genus + family b2 = 24%, SE = 9%) (Figure S3, Table S1.3). These results imply that these datasets are essentially orthogonal, as the b2 from the multiple ORM analysis approximated the sum of the individual ORM analyses (Figure S3).Taxonomic and functional datasets. We additionally explored whether combinations of various taxonomic and functional ORMs could improve the variance explained in traits. Combining taxonomic ORMs increased b2 (e.g., dietary PC1 – common species b2 = 13%, SE = 8%; common genera b2 = 14%, SE = 8%; common families b2 = 4%, SE = 4%, all of these combined multiple ORM b2 = 24%, SE = 9%), as did combinations of functional ORMs (e.g., for rBSC – common EClevel4 b2 = 34%, SE = 12%; common TCDB b2 = 3%, SE = 5%; common MetaCyc pathways b2 = 16%, SE = 9%, all of these combined multiple ORM b2 = 35%, SE = 13%), as well as combinations of taxonomic and functional ORMs (e.g., for age – common species b2 = 33%, SE = 7%; rare species b2 = 53%, SE = 9%; EClevel4 b2 = 56%, SE = 10%; TCDB b2 = 42%, SE = 10%; MetaCyc pathway b2 = 40%, SE = 9%; all of these combined multiple ORM b2 = 79%, SE = 10%) (Figure S3, Table S1.3).Informativeness of functional datasets. We found that some functional datasets captured similar b2 with fewer features: e.g, for dietary PC1 b2 = 48% (SE = 13%) based on 1,834 common features (combining Enzyme Commission Level 4 and TCDB variables) versus b2 = 48% (SE = 15%, p = 8.1e-5) based on n = 251,617 common genes.
Prediction using non-additive models
We note that estimation of b2 as the upper limit of predictability assumes an additive model. Hence, we also estimated variance explained by an adaboost (Freund and Schapire, 1997) model, which does not assume additivity.Using adaboost, we performed 5-fold cross-validation with 1000 iterations to determine the accuracy of predicting ASD diagnosis based on n = 607 bacteria (clr-transformed counts). We restricted this dataset to n = 99 ASD and n = 96 UNR participants to avoid confounding of relatedness through random sampling of sibling pairs. To generate our training dataset, we sampled the same number of ASD and UNR children (n = 76) across each iteration. For adaboost, we used the JOUSBoost implementation in R, and used the settings maxdepth = 10, and n_rounds = 100. We found that diagnosis was, on average, correctly classified 53% (SD = 7%) of the time (Data S2). This suggests that the microbiome data has negligible ability to predict ASD diagnosis, when assuming a non-additive model.To benchmark this analysis, we also repeated this same 1000-iteration 5-fold cross-validation model using age as the dependent variable, given the b2 results indicate this phenotype should be predicted by the microbiome. In this analysis, we again restricted to n = 99 ASD and n = 96 UNR participants, and input n = 607 bacteria as predictive features. To generate this training dataset, we randomly sampled 4/5 of the dataset for each iteration. We used the same JOUSBoost::adaboost settings. For simplicity of evaluating predictive ability, we calculated accuracy as the model’s ability to discriminate age as younger or older than 10 (ie. a binary variable). We found that age (younger or older than 10) was on average correctly classified 62% (SD = 7%) of the time, indicating age was more predictable than ASD diagnosis, as expected based on the OREML analyses (Data S2).
Differential abundance analysis
Species-level taxonomic data
We used ANCOMv2.1 (Mandal et al., 2015) (implemented in R: https://github.com/FrederickHuangLin/ANCOM) for differential abundance analysis as it is robust to statistical assumptions (Weiss et al., 2017). Briefly, for each focal feature (i.e., species or gene), ANCOMv2.1 performs additive-log-ratio transformation (offset = 1), runs statistical tests for each transformed taxa against selected predictors, then counts how many pairwise tests pass the Benjamini-Hochberg false discovery rate threshold (0.05), to generate a W-statistic for that initial feature (representing the number of null hypotheses rejected). Then, taxa are declared differentially abundant based on the W-statistic quantile, where a detection threshold >0.7 is commonly used to declare significant differential abundance. In this analysis, we used zero outlier cut-off, as we had already excluded variables with > 10 non-zero values (zero_cut argument), and did not apply a library cut-off (lib_cut) as all samples had over 4.8 million reads. Overall, we tested 607 species.We performed extensive sensitivity analyses for species-level differential abundance: without covariates, restricting to AAB participants with age ≤10, using a different taxonomic and functional profiling pipeline (MetaPhlAn2 (Truong et al., 2015), excluding participants with current antibiotic use, removing the SIB group, considering only sibling pairs, and excluding participants with fewer then 7M reads after rarefaction (Methods S1).
Gene-level functional data
We leveraged the species-level results to narrow our ANCOM analysis of differentially abundant genes to the set that are encoded by differentially abundant species. This was possible as we had access to per-individual tables detailing the abundance of each gene within each species. We used the same filters as for the species-level taxonomic data, ultimately testing 4,950 genes.
Romboutsia timonensis gene annotations
Here, we provide more detail on the annotations for the differentially-abundant Romboutsia timonensis-encoded genes, all of which had significantly lower abundance in the ASD group. All 6 significant (detection threshold > 0.7) genes had UniRef90 identifiers, all of which mapped to Romboutisa genera in the UniProtKB database, consistent with our focus on Romboutsia timonensis (Table S2.12). The associated proteins were: Aspartate-semialdehyde dehydrogenase (ASA dehydrogenase) (ASADH) (EC 1.2.1.11) (Aspartate-beta-semialdehyde dehydrogenase) – amino acid biosynthesis (L-lysine, L-methionine, L-threonine), which is involved in amino acid biosynthesis; Amidophosphoribosyltransferase (ATase) (EC 2.4.2.14) (Glutamine phosphoribosylpyrophosphate amidotransferase) (GPATase)– involved in purine metabolism and metabolism of L-glutamate to L-glutamine; Thymidylate kinase – involved in the DNA synthesis (and particularly pyrimidine) pathway; Galactokinase (EC 2.7.1.6) (Galactose kinase) – involved in galactose metabolism; Germination protease (EC 3.4.24.78) (GPR endopeptidase) (Germination proteinase) (Spore protease) – involved in bacterial spore germination; and Ribonuclease 3 (EC 3.1.26.3) (Ribonuclease III) (RNase III) – involved in dsDNA digestion.
Differential-abundance analyses for IQ-DQ
We tested for association of species abundance with IQ-DQ composite score, finding poor concordance between analyses when the SIB group was included versus excluded (with covariates: age, sex, dietary PC1-3). Only Bifidobacterium sp002742445 exceeded detection threshold > 0.7 in the ASD versus SIB+UNR analysis, and was fifth-most differentially abundant in the ASD versus UNR analysis (Figure 3E-3F, Table S2.13-S2.14). We note the limitations of this analysis: ASD diagnostic status being associated with lower IQ, and the inherent limitations in applying IQ measures (which rely on verbal ability) to people diagnosed with ASD.
Ordination
We performed ordination of the taxonomic data using the mixOmics package (Rohart et al., 2017). This involved principal component analysis after removing variables with < 10 non-zero values and performing clr-transformation.
Diversity analysis
We calculated alpha-diversity using the Shannon Index, to account for both richness and evenness in the dataset. As a sensitivity analysis, we also tested the effect of using richness and Simpson Index to measure alpha-diversity (Methods S1, Table S3). To calculate beta-diversity, we generated a weighted Unifrac index matrix, calculated using a phylogenetic tree for 1,731 bacterial species (i.e., not including archaea). We used this beta-diversity matrix to quantify differences in microbiome profiles between ASD, SIB and UNR groups using PERMANOVA (Anderson, 2008) (adjusted for age and sex) and PERMDISP2 (Anderson et al., 2006) (noting that PERMDISP2 does not enable adjustment for covariates).
Taxonomic diversity and ASD diagnosis
There were negligible group differences in alpha-diversity at the species-level using multiple alpha-diversity metrics (statistics for ANOVA tests: Shannon index p = 0.13, richness p = 0.13, Simpson index p = 0.08), including after adjusting for age, sex and dietary PC1-3 (statistics for ANOVA tests: Shannon index p = 0.31, richness p = 0.34, Simpson index p = 0.37) (Figure S5B), and in extensive sensitivity analyses (Methods S1, Table S3). Regressing species diversity against age, sex and dietary PC1-3, higher alpha-diversity was associated with older age (b = 0.037, p = 1.3e-7) and dietary PC2, representing a high-dairy diet (b = 0.042, p = 0.022). There were no group differences in beta-diversity (weighted Unifrac) (PERMANOVA p = 0.20; Figure S5A), nor difference in dispersion in ASD compared to SIB and UNR groups (PERMDISP2 p = 0.85; Figure S5A).
Diversity analysis with added covariates
We tested whether the relationship between dietary and taxonomic diversity changed when including dietary PCs and rBSC as covariates (n = 233 due to data missingness). Although dietary PC2-3 were significantly associated with dietary diversity and improved the model R2, the dietary-taxonomic diversity relationship was still marginally significant (Figure S7A-S7B).We also tested the effect of including energy intake as a covariate, and found it had no significant effect on the relationship between dietary and taxonomic diversity (Figure S7C).
Clinical and PGS relationships with diversity
We used linear models to test for associations between diversity measures (dietary and taxonomic) and potentially-related phenotypic measures (rBSC, ADOS2/G comparison score, ADOS-2/G repetitive and restricted behavior (RRB) Calibrated Severity Score, ADOS-2/G social affect Calibrated Severity Score, Social Responsiveness Scale t-score (SRS), Short Sensory Profile raw sensor score (SSP)), and biological measures (polygenic scores (PGS) for ASD (Grove et al., 2019), ADHD-ASD-TS cross-trait (Yang et al., 2021) and neuroticism (Nagel et al., 2018), and CD4+ T cell proportions). To account for multiple-testing, we performed Benjamini-Hochberg FDR correction (Table S3). We performed extensive sensitivity analyses to test these models (Methods S1). Our multiple-testing adjustment strategy was based on the number of hypotheses we aimed to test; these hypotheses are explicitly enumerated in Table S3, and sensitivity analyses for these hypotheses are also explained therein.
Other potential upstream mediators of diversity
Neuroticism. Others have proposed links between anxiety, neuroticism, and the microbiota (Kim et al., 2018; Yang et al., 2019), and there is evidence for genetic correlation between functional gastrointestinal disorders (such as irritable bowel syndrome) and both anxiety and ASD (nominally-significant) (Wu et al., 2021). Furthermore, both anxiety (Simonoff et al., 2008; Sukhodolsky et al., 2008) and gastrointestinal issues (Chaidez et al., 2014; McElhanon et al., 2014) often co-occur with ASD. In the absence of comprehensive anxiety trait measures in the AAB and QTAB, because the anxiety GWAS is underpowered, and as there is genetic correlation between ASD and neuroticism (Grove et al., 2019) we tested for association between neuroticism PGS and both dietary and taxonomic diversity. We found no association with either dietary (p = 0.94) or taxonomic (p = 0.43) diversity (Figure 6K–6L). Nonetheless, further investigation into the role of anxiety is warranted.CD4+ T cell proportions. We additionally assessed associations between immune cell proportions and dietary or taxonomic diversity, focusing on CD4+ T cell proportions as this cell type has been associated with atopy, and appears to have important roles in the gastrointestinal immune system through interactions with the microbiota. The rationale for this analysis is that increased rates of food allergies and intolerances have been reported in autism (Bresnahan et al., 2015), which may affect dietary diversity via both biological and social (e.g., parental intervention) mechanisms. Biologically, food allergies and intolerances may cause gastrointestinal symptoms such as abdominal pain, diarrhea, vomiting, and feeding issues. Socially, parents who are aware of this possibility may intentionally restrict their children’s diet to identify whether this improves behavior. Using data on clr-transformed CD4+ T cell proportions in n = 150 individuals, we found no evidence for associations with either dietary (b = 1.3e-2, p = 0.44) or taxonomic (t = −0.20, p = 0.10) diversity (Figure 6I and 6J), acknowledging that power is likely low.
Estimation of metabolite production potential
We estimated short-chain fatty acid metabolite production potential using the functional MetaCyc pathway dataset. For a given metabolite, we identified contributing pathways (Table S5) and summed their read counts.
Plots
We used the R packages ggplot2 (Wickham, 2016) and ggstatsplot (Patil, 2021) for plotting.
Data and code availability
The AAB datasets supporting the conclusions of this article are available by application to the Australian Autism Biobank within the Cooperative Research Centre for Living with Autism (Autism CRC): https://www.autismcrc.com.au/biobank. The QTAB dataset used in these analyses is available with mediated access: UQ eSpace: https://espace.library.uq.edu.au/view/UQ:e803a68 . Code is publicly available at https://zenodo.org/record/5558047. Any additional information required to reanalyse the data reported in this paper is available from the Lead Contact upon request.
Acknowledgments
We thank all participants and families in the Australian Autism Biobank, without whom this study would not have been possible. The data used in this project were provided by the Co-operative Research Centre for Living with Autism (Autism CRC) with appropriate ethics approval. The Autism CRC is established and supported under the Australian Government’s Cooperative Research Centre Program. We thank Dr. Natalie Silove and the Child Development Unit team at Sydney Children’s Hospital, Westmead; the staff at KU Marcia Burgess Autism Specific Early Learning and Care center and KU Children’s Services; Dr. Nancy Sadka from the Olga Tennison Autism Research Centre at La Trobe University for her integral role in recruitment at the Victorian site; and Dr. Anja Ravine and the pathology services at Pathwest for their contributions in collecting and processing blood samples in Perth. We thank Felicity Rose for assistance with the AAB. The QTAB Project thanks the twins and their families for their willingness to contribute to research and for giving their time so generously. Special thanks to the many research assistants in the QTAB team, and to the Human Studies Unit, Program in Complex Trait Genetics, IMB, UQ, for the processing of QTAB samples and the Queensland Twin Registry Study. We thank Zhiyu Yang for providing GWAS summary statistics for the ADHD-ASD-TS PGS analysis and Kylie Ellis for assisting with metagenomics analysis. The authors acknowledge the financial support of the Autism CRC. QTAB was facilitated through access to Twins Research Australia, a national resource supported by a Centre of Research Excellence Grant ( 1078102 ) from the Australian National Health and Medical Research Council . We also acknowledge funding support from the Australian National Health and Medical Research Council ( 1103418 and 1127440 to J.G., 1078901 and 1173790 to N.R.W., 1113400 to N.R.W. and P.M.V., 1077966 and 1173896 to A.J.O.W., and 10787561 to M.J.W.), the Australian Research Council ( FT200100837 to A.F.M. and FL180100072 to P.M.V.), and the University of Queensland (RTP Stipend and Tuition Fee Offset as well as the Sam and Marion Frazer HDR Top-up Scholarship in Neurological Disease to C.X.Y.). This work was supported by Mater Research and the Mater Foundation and was carried out in part at the Translational Research Institute (TRI), which is supported by a grant from the Australian Government .
Author contributions
The Australian Autism Biobank was conceived and designed by A.J.O.W., C.D., V.E., H.S.H., G.A.A., P.A.D., J.G., R.G., A.K.H., L.P.L,. A.M., L.W., and N.R.W. Australian Autism Biobank patient recruitment, assessments, and data collection were led by A.J.O.W., V.E., C.D., M.L.F., H.S.H., and P.A.D. G.A.A., D.C., R.G., C.H., A.H., H.H., R.J., F.K., L.P.L., J.L., M.L.F., A.M., N.E.M., M.M., and M.N. contributed to the coordination, data collection, and management through the four sites. Samples were processed by L.W. and T.M., overseen by A.K.H. The Queensland Twin Adolescent Brain Project was conceived and designed by M.J.W., G.I.d.Z., P.M.T., and K.L.M. with N.K.H. and L.T.S. contributing to ongoing project design and management, A.K.H. overseeing the management of biological samples, and J.L.M. and L.N. processing samples. J.G., C.X.Y., and A.K.H. designed the metagenomics study. C.X.Y. performed analyses with metagenomics QC and profiling by D.L.A.W. and L.K. and with advice from R.R., G.A.A., A.F.M., G.W.T., and G.H. J.G. and N.R.W. supervised the analyses. C.X.Y., J.G., and N.R.W. wrote the manuscript with critical input from all authors.
Declaration of interests
David L.A. Wood and Lutz Krause are employees of Microba Life Sciences. Gene W. Tyson is a co-founder and director of Microba Life Sciences. Gerald Holtmann is on the advisory board of Servatus Biopharmaceuticals. The other authors declare no competing interests.
Inclusion and diversity
We worked to ensure that the study questionnaires were prepared in an inclusive way. We worked to ensure ethnic or other types of diversity in the recruitment of human subjects. We worked to ensure gender balance in the recruitment of human subjects. One or more of the authors of this paper self-identifies as an underrepresented ethnic minority in science. One or more of the authors of this paper received support from a program designed to increase minority representation in science.