Comparison of Fecal Microbiota in Children with Autism Spectrum Disorders and Neurotypical Siblings in the Simons Simplex Collection

Families were recruited via a registry called the Simons Simplex Community through the Interactive Autism Network (SSC@IAN), which is composed of simplex families originally recruited to the Simons Simplex Collection that were willing to be contacted through the Interactive Autism Network for additional studies. The SSC@IAN provided us with contact information for 245 potentially eligible simplex families that expressed their willingness to be contacted about this study. Of the 245 potentially eligible SSC families with 339 children, 107 families with 145 children consented to participate with parental written consents and written child assents for all participating children, and 66 families with 103 children (66 ASD probands and 37 NT siblings) completed the study. Exclusion criteria for ASD probands and NT siblings in this study were age <7 or age >14 years, a past or current history of GI diseases or other serious medical problems, and use of antibiotics within one month of collection of the stool sample. The procedures of the present study were reviewed and approved by the Institutional Review Board of Stony Brook University (IRB 313466)

Families who consented to participate in the study were asked to complete the Pediatric ROME III Version (QPGS-RIII), the GI symptoms form in the Simon Simplex Collection Medical History (v 4.0) [31], the Child Behavior Checklist (CBCL/6-18) [33–37], a one-week food diary for the ASD proband and/or the NT sibling prior to collecting a fecal sample, and provide current information on the child’s height and weight. The families were asked to submit the stool samples after the subjects were off probiotics and antibiotics for at least one month.

Functional GI disorders (FGID) were assessed using the ROME III process [25]. Information on GI symptoms in both ASD and NT siblings was also collected using the GI symptoms section of the SSC Medical History form [32]. The GI symptoms section of the SSC Medical History form included the seven symptom categories: constipation, diarrhea, abnormal stool smell, bloating, excessive gas, severe abdominal pain, and vomiting. The parents were asked to rank whether these categories were problematic as “yes”, “not sure” and “no”. Many of these symptom categories in the SSC form overlap with the following symptoms included in the calculation of the GI severity index score [27]: constipation, diarrhea, stool smell, flatulence and abdominal pain, but did not include “unexplained daytime irritability”, “nighttime awakening” or “abdominal tenderness during exam”.

The CBCL/6-18 is a parent report form to screen for emotional, behavioral, and social problems in children ages 6–18 [33]. This instrument as well as the earlier version of this instrument has been previously used by several research groups studying school-aged children with ASD [34–37]. The data from ASD and NT siblings was entered and analyzed using the AGES 6–18 ASEBA-PC MODULE (ASEBA) to yield norm-referenced T-scores (M = 50, SD = 10). The CBCL contains eight syndrome scales: Anxious/Depressed, Withdrawn/Depressed, Somatic Complaints, Rule Breaking Behavior, Aggressive Behavior. The Internalizing Domain score is a measure of emotional problems and contains three syndrome scales: Anxious/Depressed, Withdrawn/Depressed, and Somatic Complaints. The Externalizing Domain score is a measure of behavioral problems and contains the Rule Breaking Behavior and Aggressive Behavior syndrome scales.

The parents were asked to record the daily dietary intake of ASD and NT siblings on the Stony Brook University Medical Center Department of Pediatrics Food Diary/Calorie Count sheet for 7 days prior to the collection of the stool sample. The data were entered and analyzed using Nutritionist Pro nutrient analysis software version 5.2.0 (Axxya Systems, Nutritionist Pro, Stafford, TX). We focused on analyzing the following six macronutrient categories: total kcal/day, total protein g/day, total carbohydrates g/day, total fat g/day, total dietary fiber/day, total sugar g /day.

Fecal samples were collected in a stool “hat” placed in the toilet and following evacuation 2 ml of stool were immediately transferred, using the spoon built into the lid into a Para-Pak Clean vial (Meridian Bioscience, Cincinnati, OH) that contained 10 ml of RNAlater for metagenomics studies and a Para-Pak Clean vial with no preservative. The fecal samples were shipped in cold packs overnight to our laboratory. Fecal DNA was extracted from stool samples immediately upon arrival using ZR Fecal DNA MiniPre (Zymo Research Corporation, Irvine, CA) and following manufacturer’s protocol. Fecal DNA and stool samples were stored in -80°C.

QPCR assays were performed using established primers for total bacteria (forward, 5’- GTG STG CAY GGY TGT CGT CA-3’ and reverse 5’- ACG TCR TCC MCA CCT TCC TC-3’) [38], Sutterella subgroup (forward, 5’-CGC GAA AAA CCT TAC CTA GCC-3’ and reverse 5’- GAC GTG TGA GGC CCT AGC C-3’) [10], Bacteroidetes subgroup (forward 5'-AACGCTAGCTACAGGCTT-3’ and reverse 5'-CCAATGTGGGGGACCTTC-3’) [39]; Prevotella (forward, 5’-CACCAAGGCGACGATCA-3’ and reverse 5’-GGATAACGCYGGACCT-3’) [40], C. coccoides-E. rectales subgroup (forward, 5’—CGGTACCTGACTAAGAAGC-3’and reverse 5’-AGTTT(C/T)ATTCTTGCGAACG-3’) [41], Faecalibacterium prausnitzii (forward 5’-CCC TTC AGT GCC GCA GT-3’ and reverse 5’-GTC GCA GGA TGT CAA GAC-3’) [42], Escherichia coli subgroup (forward, 5’- GTT AAT ACC TTT GCT CAT TGA-3’ and reverse, 5’- ACC AGG GTA TCT AAT CCT GTT-3’) [42]. The log₂ transformation of the relative abundance of each bacterial subgroup was measured by ΔCt = Ct (threshold cycle) _{total bacteria}−Ct _subgroup. All assays were carried out in triplicate. Plasmid quantification standards were prepared from representative clones of the target organisms to insure that the assays were conducted within the linear range and with similar slopes.

Fecal bacterial profiles were determined using broad-range amplification of the ~300 base pair (bp) V1V2 variable region using primers 27FYM (AGAGTTTGATYMTGGCTCAG) and 338R (TGCTGCCTCCCGTAGGAGT) and the ~500 bp V1V3 region using primers 27FYM (AGAGTTTGATYMTGGCTCAG) and 534R (ATTACCGCGGCKGCTGG) [43,44]. Sequencing was performed on the Illumina MiSeq platform using 600-cycle version 3 reagent kits. Demultiplexed, paired end sequence data were deposited in the NCBI Short Read Archive under BioProject Accession Number: PRJNA282013 (www.ncbi.nlm.nih.gov/bioproject/PRJNA2282203↗). The sorted paired reads were assembled using phrap [45] and paired reads that did not assemble were discarded. Assembled sequence ends were trimmed over a moving window of 5 nucleotides until average quality met or exceeded 20. Trimmed contigs with more than 1 ambiguity or shorter than 200 nucleotides were discarded. Potential chimeras identified with Uchime (usearch6.0.203_i86linux32) [46] using the Schloss [47] Silva reference sequences were removed from subsequent analyses. Assembled sequences were aligned and classified with SINA (1.2.11) using the 418,497 bacterial sequences in Silva 115NR99 as reference configured to yield the Silva taxonomy [48,49]. Operational taxonomic units (OTUs) were produced by clustering sequences with identical taxonomic assignments. OTU counts were normalized between samples by dividing sequence counts by the total number of sequences generated per sample to calculate the relative abundance. Phylum-level OTU tables were generated by collapsing lower level OTUs into higher-level categories. OTUs with a maximal relative abundance <0.0001 and with a prevalence <0.01 were culled. The remaining OTU relative abundances were then transformed using the square root function.

Comparisons of gender, FGIDs, and SSC GI symptoms between ASD and NT siblings were conducted using the unconditional Fisher exact test with one tail [50]. The threshold significance was set at P < 0.05.

The ASD and NT siblings were further subdivided with respect to whether they also carried a diagnosis of an FGID based on the ROME III process. Study characteristics including body mass index (BMI), CBCL scores, dietary scores and the qPCR ΔCt scores were compared between the following four groups: 1) ASD with FGID, 2) ASD without FGID, 3) NT siblings with FGID and 4) NT siblings without FGID. Outcome variables including CBCL, dietary macronutrients, qPCR ΔCt were also modeled using mixed effect models with ASD, FGID and ASD*FGID as the three fixed effects, and family as random effects. Each continuous variable was fitted to a linear model while each categorical variable was fitted to a logistic model. Comparison of continuous variables between family matched ASD and NT siblings was conducted using two-tailed Student’s T test or the Wilcoxon signed rank test with the threshold of significance set as P < 0.05.

Alpha diversity indices (e.g. Chao1, Shannon complexity H, Shannon Evenness H/Ho) were calculated using Explicet [51] for each of the four groups with 10,000 replicate resamplings: 1) ASD with FGID, 2) ASD without FGID, 3) NT siblings with FGID and 4) NT siblings without FGID.

Beta diversity was compared between ASD and NT siblings using the adonis function in the R vegan package as previously described [52,53] at the phyla, family and genus level. This function uses a non-parametric multivariate analysis of variance test (PERMANOVA) [54], and was applied using Bray-Curtis, Jaccard and Morisita-Horn indices as distance measurements [55]. Family (i.e. whether ASD or NT siblings belonged to the same family) was selected as a strata variable to restrict permutations as previously described [56]. Overall microbial composition binned as phyla, was also compared between family matched ASD and NT siblings using the permutation Hotelling T test [57].

Paired two class comparisons of the transformed (square root) relative abundances of individual OTUs (phyla and genera level) between ASD probands and NT siblings from the same family were conducted using the Wilcoxon signed rank test with the threshold set as FDR < 0.05.

Because over-dispersion is often observed in microbiome sequence count data [58], a negative binomial regression model was also used to assess the effect of the ASD phenotype, FGID phenotype and first order interactions between ASD and FGID (ASD*FGID) on each individual OTU, in both paired and unpaired ASD and NT samples as follows: log(μijk)=(βi0)k+β1kASDij+β2kFGIDij+β3kASD*FGIDij+log(Yij.)(βi0)k=b0k+bikI{familyij=i}Yijk∼NB(μijk,ϕk). where Y_ijk denotes the raw sequence count of OTU k for child j in family i. In this model we treat family, which may reflect genetic and possibly dietary factors, as a zero mean random coefficient. Coefficients b₀, β₁, β₂ and β₃ are fixed constants representing grand mean, the effect of ASD, FGID and the interaction of ASD and FGID respectively. The log total sequence count (Yij.=∑kYijk) for each individual is considered as an offset. Because the distribution of the raw sequence count is often skewed by a large proportion of zero counts, a zero-inflated version of the negative binomial model was also fitted to each individual OTU: The final model was selected based on the Akaike Information Criterion (AIC). The p-values for each OTU were adjusted for multiple comparisons to calculate a false discovery rate (FDR) using the Benjamini-Hochberg procedure [59]. The threshold of significance for ASD, FGID or ASD*FGID first order interactions was set at FDR < 0.1. The logistic and linear mixed models were built with R package glmmADMB [60] and lme4-R [61] respectively.

Comparative analysis of the parallel square root transformed relative abundance [62] measurements of Bacteroidetes phylum, the Sutterella genus and the Prevotella genus by MiSeq V1V2, MiSeq V1V3 and targeted qPCR was conducted using latent structural equation modeling [63]. Because of the large proportion of zero counts for Sutterella and Provetella genera, samples with zero counts were excluded from the analysis. To avoid the complication caused by correlated samples, SEM models were fitted with grouping variable ASD and the constraint that loadings are equal for different groups. The corresponding model for each group is: The corresponding model for each group is: qPCRij=a1+b1ξijk+δ1V1V2ij=a2+b2ξijk+δ2V1V3ij=a3+b3ξijk+δ3 Here ξ_ij is the underlying true bacteria relative abundance in the square root scale, while δ₁, δ₂ and δ₃ are the corresponding measurement errors. In order to render the model identifiable, we chose qPCR as the reference level with a₁ = 0 and b₁ = 1. The model was fitted with the R package lavaan [64].

A total of 59 ASD children and 44 NT siblings from 66 families completed the study (Table 1). Matched ASD and NT sibling pairs were from 37 families. In some instances NT siblings of ASD probands were recruited, but not the ASD proband s because their age exceeded 14 years. The ASD children were predominantly male. The gender distribution for the NT siblings was roughly equal. The mean age in both ASD children and NT siblings was 10 years, and they were predominantly of Caucasian descent. There are three major groups of FGIDs (see Table 1) [24]: H1) Vomiting and aerophagia, H2). Abdominal-pain related FGIDs, and H3) Constipation and fecal incontinence. Group H1 is further subdivided into three FGIDs: H1a) Adolescent rumination syndrome; H1b) Cyclic vomiting syndrome; H1c) Aerophagia. Group H2 is further subdivided into four FGIDS: H2a) Functional dyspepsia, H2b) Irritable bowel syndrome (IBS), H2c) abdominal migraine, H2d) Childhood functional abdominal pain with a subset termed H2d1. Childhood functional abdominal pain syndrome. Subjects classified within the H2d1 subset have additional somatic symptoms such as headache, limb pain or difficulty sleeping.

Twenty-five of 59 (42%) ASD children and 13 of 44 (30%) NT siblings were diagnosed with at least one FGID. Fifteen ASD children and 13 NT siblings in 37 families met ROME III criteria for at least one FGID, with concordance for FGIDs in 6 (18%) ASD/NT sibling pairs. Some ASD and NT siblings were diagnosed with more than one FGID. ASD children exhibited a higher prevalence of functional constipation in compared to NT siblings (P = 0.035) when both matched and unmatched individuals were included in the analysis. We detected no significant difference in the prevalence of children with at least one FGID in family matched ASD/NT sibling pairs. We also detected no significant difference in the prevalence of specific FGIDs such as functional constipation in the family matched ASD/NT sibling pairs.

When the SSC registry was established, medical history forms including one on GI symptoms were completed by the parents on the ASD children but not the NT siblings. The same SSC GI symptoms form was then completed (~ 2–3 years later) by the parents for the ASD and NT children for the current study. As shown in Table 2, there was no significant change in the prevalence of the GI symptoms occurred over time. Furthermore, the SSC results indicate that ASD children exhibited a higher prevalence of constipation compared to NT siblings (P = 0.05).

As shown in Table 3and S1 Table, the ASD severity score, ADOS-CSS, did not differ significantly between the ASD w FGID group and the ASD w/o FGID group. Mean total, internalizing, externalizing and subscale CBCL scores were significantly increased (p < 0.05) in the ASD children in both the ASD/NT sibling pairs and in comparison of all the ASD with NT children (matched and unmatched). These CBCL scores were also significantly increased (p < 0.05) in NT siblings with FGID compared to NT siblings without FGID (Table 3). There was no significant difference in the mean CBCL scores between ASD children with FGID compared with ASD children without FGID.

No significant differences (p < 0.05) between the four study groups with respect to daily intake of macronutrients (calories, protein, fats, carbohydrates, sugars or dietary fiber) were noted (Table 3). Only 4 children were reported to adhere to a gluten and/or casein free diet. The study participants had BMIs that were within the normal range and were not obese.

Targeted qPCR assays of fecal samples were conducted for Sutterella, Prevotella and total Bacteroidetes subgroups, because of previous reports that the relative abundance of these subgroups differed between ASD and NT children [7,8,10–14]. QPCR assays targeting the C. coccoides-E-rectales group, Faecalibacterium prausnitzii and Escherichia coli were also performed because altered relative abundances have been reported previously in GI disorders such as inflammatory bowel diseases [23]. The mean ΔCt values of these subgroups, which represent a log₂ transformation of the relative abundance, for 1) ASD w FGID, 2) ASD w/o FGID, 3) NT w FGID and 4) NT w/o FGID are shown in Table 4. No significant difference was observed for any of these bacterial subgroups in the 37 ASD/NT sibling pairs (see Methods). Furthermore, in mixed models, no significant effect of ASD, FGID or ASD*FGID was observed in an analysis of all 59 ASD and 44 NT siblings.

We performed Illumina sequencing of the V1V3 hypervariable region of the bacterial 16S rRNA gene to facilitate comparisons of this data with previous 16S rRNA 454 pyrosequencing datasets of the V1V3 hypervariable region [16,23]. A total of 11,050,770 high-quality sequences were generated (median 1.1E+05 sequences/sample, IQR 8.0e+04 to 1.2E+05). All libraries had a Good’s coverage score ≥ 99.9% at the rarefaction point of 15,000 sequences, indicating deep sequence coverage of the intestinal microbiome was achieved for each sample. To corroborate results generated through V1V3 sequencing, we also performed Illumina sequencing of the V1V2 hypervariable region. A total of 14,783,994 high-quality sequences were generated (median 1.4E+05 sequences/sample, IQR 1.3e+05 to 1.6E+05; minimum Good’s coverage = 99.9%). No significant effect of ASD, FGID or ASD*FGID was observed on OTU complexity or richness, two measures of α-diversity, Shannon diversity and Chao 1 (see Table 5for mean values). No significant difference (p < 0.05) was noted in β-diversity between ASD and NT sibling bacterial communities analyzed at the phyla, family or genus level and using Bray-Curtis, Jaccard or Morisita-Horn indices as dissimilarity measurements.

The relative abundances of the major four phyla, (Firmicutes, Bacteroidetes, Actinobacteria, Proteobacteria) for the ASD children with FGID, ASD children without FGID, NT siblings with FGID and NT siblings without FGID groups for the V1V2 and V1V3 datasets are shown in Fig 1(see S2 Table). Phyla of lower abundance, including RF3, Tenericutes, Cyanobacteria, Verrucomicrobia, Lentisphaerae, and Fusobacterium were also detected (see S3 Table). Permutational Hotelling T test revealed no significant difference in the overall microbial composition when binned at the phyla level, between ASD children and NT siblings. With the exception of the Verrucomicrobia phylum, the relative abundances of the phyla were similar between the V1V2 and V1V3 data sets. Exploratory analyses of individual phyla detected no significant effect (FDR <0.1) of ASD, FGID or ASD*FGID first order interactions for any of the individual phyla, including Bacteroidetes, using a negative binomial model. At the genus level, 359 operational taxonomic units (OTU) were detected for the V1V2 dataset and 189 OTUs for the V1-V3 datasets, respectively. No significant associations were detected between the relative abundance of the Sutterella or Prevotella genera and ASD, FGID or ASD*FGID first order interactions (see Table 6) in either the V1V2 or V1V3 datasets. Exploratory analyses of individual genera revealed a significant effect (FDR < 0.1) of ASD and ASD*FGID first order interactions on the Cyanobacteria/Chloroplast genus for both the V1V2 and/or V1V3 datasets using a zero-inflated negative binomial model (see Table 6). The mean relative abundance of this genus was increased in the ASD-with FGID group compared to those of the other three categories (ASD without FGID, NT siblings with FGID and NT siblings without FGID, see S4 Table for mean relative abundances). 16S RNA sequence analyses of the human gut microbiota have previously detected sequences binned within the Cyanobacteria phylum. These sequences could represent chloroplast 16S rRNA genomic sequences present in dietary plant foods or novel uncultivated bacterial species. Metagenomic studies have detected a novel nonphotosynthetic taxa termed “Melainabacteria” present in human gut flora that represents either a sibling phylum or a class within the Cyanobacteria phylum [65,66]. The Silva database used in the current study, however, binned the Melainabacteria sequences in a taxon that was distinct from the Chloroplast taxon. Inspection of the daily food diaries of all the subjects revealed that two ASD children with FGID (functional constipation), who had high relative abundances of this genus and consequently dominated the signal associated with ASD with FGID category, were fed chia seeds added to smoothies. None of the other subjects consumed chia seeds except one ASD child who ingested a serving of mixed grain bread that included chia grain. Of note, chia seeds have been touted as an autism cure on the internet [67,68].

For the V1V2 dataset only, significant associations were observed for 1) ASD and ASD*FGID first order interactions with Firmicutes/Asteroleplasma, 2) ASD, FGID and ASD*FGID first order interactions with Proteobacteria/Thalassospira, and 3) FGID and ASD*FGID first order interactions with Proteobacteria/Burkholderia (see Table 5and S4 Table). The Firmicutes/Asteroleplasma genus appears to represent a group of uncultured wall-less mycobacterium like organisms within the class Erysipelotrichia [69]. The Proteobacteria/Thalassospira genus appears to represent alpha-proteobacteria present in bacterioplankton and a number of uncultivated organisms [70]. The sequences binned in this Proteobacteria/Burkholderia taxon represent unclassified sequences within the Burkholderia order that did not bin within the Comamonadaceae family (see below).

For the V1V3 dataset only, significant associations were observed for 1) ASD and ASD*FGID first order interactions with the Proteobacteria/Comamonadaceae, 2) ASD, FGID and ASD*FGID first order interactions with Fusobacteria.Fusobacteriales, 3) ASD*FGID first order interactions with Bacteroidetes/Prevotellaceae, and 4) ASD and ASD*FGID interactions with Actinobacteria/Mobiluncus. These associations, however, did not reach significance in the V1V2 dataset. Proteobacteria/Comamonadaceae is a large and diverse bacterial family within the Burkholderiales order (see above). The sequences binned within this Fusobacteria/Fusobacteriales taxon represent unclassified sequences within the Fusobacteriales family that did not bin within the Fusobacterium genus. The sequences binned within the Bacteroidetes/Prevotellaceae taxon represent unclassified sequences within the Prevotellaceae family that did not bin within the Prevotella genus. The Actinobacteria/Mobiluncus genus has previously been associated with bacterial vaginosis [71].

Comparisons of the relative abundances of the Bacteroidetes phylum, and the Sutterella and Prevotella genera reveal high correlations between MiSeq V1V2 and V1V3 (ρ≈0.9) and small proportional biases (S5 Table and S1 Fig). The different reverse primers, 338R and 534R, used for amplifying respectively the V1V2 and V1V3 regions, might have subtle differences in specificity. While there were generally positively correlated measurements between the MiSeq and qPCR platforms, we observed large systematic biases particularly with the two relative low abundance genera, Sutterella and Provetella. This may reflect at least in part that the depth of sequencing was not sufficient to accurately measure relatively low abundance bacterial categories.