SNVs in ASD are related to the gut microbiome

A total of 567 exonic SNVs (Table S1) were identified in our ASD cohort using the matched TD cohort and the public databases as control (Figure 1). These SNVs shared the same mutation pattern (Figure 2(a)) and mutation signature (Figure 2(b)) with the previously reported ASD-related mutations collected in the AutismKB, a knowledgebase for the genetic evidence of autism spectrum disorder.²² We calculated the coordinates underlying variability in the host SNV genotype data using Principal Component Analysis (PCA). We then computed alpha diversity, a measure of within-sample microbial diversity (Simpson index within a sample). We found that the first coordinated of the functional SNVs (protein-altering SNVs, including non-synonymous, stop-gain, and stop-loss SNVs) was significantly correlated with alpha diversity (Figure 2(c), R = 0.45, P = .020) of the gut microbiome. The altered alpha diversity links to lots of human health disease,^23,24 including autism.¹⁵ Our results suggested a possible role for host genetic variations in shaping it. Next, we looked for the association of host genetics with the overall composition of the gut microbiome. A marginally significant correlation was observed between the first principal coordinate of functional SNVs and microbiome PC1 (Figure 2(c), R = 0.37, P = .059). However, when considering the total enriched SNVs, no significant correlation between the host genetics and microbiome diversity or composition was detected (Figure S1). These correlations suggest a potential relationship between the functional SNVs and microbiome.

Figure 1.

Schematic overview of the study

We performed whole exome sequencing of the host genome, whole genome shotgun sequencing and metabolite profile of fecal samples of 26 ASD individuals and 26 TDs. Associations were tested between genetic variations and microbes, as well as microbial genes, metabolites and pathways, and host immunity.

Figure 2.

SNVs enriched in ASD children. (A-B) The comparison of the mutation pattern (A) and mutation signature(B) between our identified SNVs and mutations collected from the AutismKb database. (C) Correlation of the first PC of host genetic data (x-axis) and alpha diversity (left panel) or first PC of the stool microbiome (right panel). (D) Circos plot of the functional SNVs enriched in ASD. The height of the points represents the number of mutant samples. Genes mapped by SNVs were labeled outside the circle. (E) Overlap of genes with the ASD enriched SNVs in our studies with those collected from other studies. (F) Functional annotation of genes with functional SNVs

The functional SNVs were mapped to 206 genes (Figure 2(d)), including 70 previously reported genes,^4,22,25,26 such as MUC4, SETD5, and ANKRD11^27-29 (Figure 2(e)). Functional enrichment analysis on these genes revealed a significant enrichment on the protein glycosylation process, innate immune response, and retinoid metabolic process (figure 2(f)). The mucins encoded by MUC family genes are O-glycosylated proteins that play an essential role in forming protective mucous barriers on epithelial surfaces,³⁰ and mucins are components of innate and adaptive immune responses to mucosal infection.³¹ Increasing evidence has shown that variants in genes encoding glycosylated extracellular proteins or enzymes (glycogene) involved in glycosylation may contribute to the etiology and pathogenesis of ASD.³² For examples, the mucin encoding gene MUC4 and the glycosyltransferase B3GALNT2.³³ Studies in mouse models of congenital disorders of glycosylation (CDGs) and behavioral phenotypes observed in CDG patients support the idea that glycogene variants either cause or contribute to the development of idiopathic ASD.³² Retinoids are vitamin A and its derivatives. The retinoid signaling plays a significant role in regulating brain functions, including neuronal differentiation, neurite growth, patterning of the anteroposterior axis of the neural tube, neurotransmitter release, and long term potentiation.³⁴ Defective retinoid signaling has been evidenced in the pathology of Alzheimer’s disease,³⁵ and Vitamin A deficiencies were reported to exacerbate symptoms in children with ASD.³⁶

SNVs correlated with bacteria taxa

Systematic screening for associations between each functional SNVs and abundance of each species was conducted using standard linear regression after adjusting for age and sex for confounding factors. There were 46 candidate associations identified (Figure 3, Table S2). Among the results, a variation in the glycosyltransferase B3GALNT2 (chr1:235659576, c.C155T, p.A52V) was associated with the increased abundance of Catenibacterium mitsuokai in a subset of ASD children. A variation of PGLYRP4 (rs148195147, c.G602A, p.R201Q) was related to Faecalibacterium prausnitzii. The PGLYRP4 is a part of the innate immune system and encodes a known antibacterial protein that recognizes peptidoglycan, a ubiquitous component of bacterial cell walls. The deficiency of PGLYRP4 leads to changes in the gut microbiota.³⁷ Other associations involved by the immunity-related genes include CD180 variation (chr5:66479217, c.C1454A, p.T485N) associated with Bacteroides dorei and Clostridium sp. CAG:7, RAET1G variation (rs4870111, c.G668A, p.R223K) associated with two Bacteroides species and Catenibacterium mitsuokai, and RTN4 variation associated with Hungatella hathewayi. CD180 is a TLR-associated molecule that is mainly expressed on macrophages, DCs, and B cells, which can also sense bacterial LPS to trigger the growth of dormant multiplemyeloma cells.³⁸ Both RAET1G and RTN4 mutations were reported in ASD cohorts.²² RAET1G encodes a member of the major histocompatibility complex (MHC) class I family of proteins. It is one of the ligands of natural killer group 2, member D (NKG2D) receptor, which functions as an activating receptor in innate and adaptive immunity.³⁹ RTN4 is necessary for immune responses triggered by nucleic acid-sensing TLRs.⁴⁰ Besides, RTN4 is also an inhibitor of axonal regeneration in the CNS.

Figure 3.

SNVs associated with microbes

SNV-microbe associations identified by regression. The SNVs were mapped to genes. The size of the node is proportional to the number of associations involved by the gene or microbe.

Also, APOB (Apolipoprotein B, rs13306187, c.G4163A, p.R1388H) and RARRES2 (Retinoic Acid Receptor Responder 2, rs147597725, c.G440A, p.S147N), two genes involved in the retinoid metabolic process, were related to the abundance of Lachnospiraceae bacterium 7_1_58FAA and Bacteroides eggerthii, respectively. Apolipoprotein B was an innate barrier against bacterial infection,⁴¹ and novel APOB mutations were observed in autism.⁴² The Bacteroides eggerthii was reported to enhance colitis in mice.⁴³ Studies have reported that the Bacteroidetes/Bacteroidales were among the key taxa related to vitamin A in ASD children.⁴⁴ Other associations include the thiamine transporters SLC19A3 (rs117864472, c.A1132G, p.I378V) and Clostridiales bacterium VE202-03, among others. The thiamine deficiency was observed in Alzheimer’s disease and contribute to synapse and neural circuit defects.⁴⁵

Previous studies have demonstrated that ASD genetic risk can act in part through expression quantitative trait loci,⁴⁶ namely, eQTL, and epigenetic regulation,⁴⁷ such as meQTL. Accordingly, we explored the association between microbes and the eQTL or meQTL SNVs. We mapped the total SNVs to the publicly available eQTL⁴⁸ and meQTL⁴⁷ maps for multiple tissues. Thirty SNVs in our ASD cohort were the eQTL or meQTL loci of neuropsychiatric disorders (Figure 4(a)). However, a limited number of these loci were associated with microbes (Figure 4(a)). An exonic variation in the non-coding RNA LOC285819 (rs34104395), which has an eQTL effect on BTN3A2 in both colon and brain tissue, was associated with Clostridium sp. L2-50. BTN3A2 is a member of the immunoglobulin superfamily, residing in the juxta-telomeric region of the major histocompatability class 1 locus. Increased expression of BTN3A2 might confer risk for schizophrenia by altering excitatory synaptic function. Previous studies have reported the influence of MHC gene variation in shaping gut microbiome composition in both mice⁴⁹ and human.⁵⁰ This observation suggested a possible genetic regulation of gut microbiota mediated by the host immune system. Besides, the eQTL loci (rs75634125) of MRPS22 (Mitochondrial Ribosomal Protein S22) was associated with the abundance of Roseburia faecis. The MRPS22 mutation was reported to be associated with the mitochondrial disorder.⁵¹

Figure 4.

The association between QTL SNVs and microbial pathways, the causal relationship between SNVs, microbes, and metabolites

(A) The association between eQTL/meQTL loci and microbes, and the genes putatively regulated by these loci. The inner labels are ASD SNVs that mapped to the eQTL or meQTL map, and the outermost labels represent genes regulated by the eQTL or meQTL loci. The blue links represent the associations between ASD SNVs and the gut microbes. (B) The associations between SNVs (gene) and microbial pathways (KEGG pathway, levele 3).

Host SNVs associated with microbial pathways and microbiome gene abundance

To assess how and to what extent host SNVs impact the functions of microbiome, associations between SNVs and microbial pathways and genes were tested. We identified 22 associations between 13 host genes and 20 microbial pathways (Figure 4b, Table S3). Several variations were associated with metabolism pathways, including variations located in MUC3A (chr7:100550525,c.C1106T, p.P369L; chr7:100550534, c.C1115G, p.S372 C) and ANKRD20A2/3 (Ankyrin Repeat Domain 20 Family Member A1, rs201420500, c.C2456A, p.S819Y) that are associated with lipid metabolism pathways, and LAMTOR1 variation (rs146341570, c.C374T, p.P125L) that associated with lysine degradation. Notably, variations in PLCG2 (rs201080992, c.C579G, p.H193Q) and ALDH1A (rs8187929, c.A529T, p.I177F) were associated with the pathway of Xenobiotics biodegradation and metabolism (caprolactam degradation and naphthalene degradation pathways). ALDH1A1 is a multifunctional enzyme with dehydrogenase, esterase, and anti-oxidant activities. Moreover, coding variation in PLGG2 implicated microglial-mediated innate immunity in Alzheimer’s disease.⁵² Toxin exposure has been epidemiologically demonstrated as one of the main etiological factors of ASD, and deficiencies in toxicant-degradation pathways were reported in ASD.⁵³ These results implicated the possible genetic factors that were underlying the microbial detoxification deficiency. Besides, a couple of gene variations were associated with the immune system. For example, variation in an ASD-related gene SSH1 (rs117900986, c.A1862G, p.N621S) was associated with the human immune disease. SSH1 encodes proteins dephosphorylate and activate the actin-binding/depolymerizing factor cofilin. It is required by NOD1 to detect bacterial-induced changes in actin dynamics leading to NF-κB activation and innate immune responses.⁵⁴ Other associations were observed between variation within NOXA1 (rs201388549, c.C526T, p.R176W), an activator of Nox2-based NADPH oxidase,⁵⁵ and the carbohydrate digestion and absorption, and variation of ANKRD42 (rs14346246, c.G1547A, p.G516D) with the wnt signaling pathway.

Furthermore, compared to the total annotated genes, the host-SNVs-associated genes were enriched in multiple metabolism pathways in the KEGG annotation (Table S4), including carbohydrate metabolism, amino acid metabolism, metabolism of cofactors, nucleotide metabolism and vitamins, and biosynthesis of other secondary metabolites. Additionally, enriched pathways include membrane transport (Figure S2A). Also, these microbial genes were enriched in cell motility for the eggNOG annotation (Figure S2B).

The causal relationship between SNV, microbiome and its metabolites

One of the critical mechanisms by which microbes regulate host physiology is through their metabolites.^15,56 To decipher how the host SNVs impact microbiome and their metabolites in ASD, a causal inference test (CIT) proposed by Millstein et al.⁵⁷ was applied to infer the causal relationship between SNV, microbes and their metabolites. Focusing on the SNVs associated with specific microbes identified in the previous section, 8 SNV-microbe-metabolite trios, comprising 4 SNV-microbe associations detected before and 8 metabolites (Figure 5(a), Table S5), were identified.

Figure 5.

The causal relationship between SNVs, microbes, and metabolites

(A) The causal relationship of SNV-microbe-metabolite trios inferred by causal inference test. The red and green links represent positive and negative assoications, respectively.(B) The abundance of Faecalibacterium prausnitzi related to PGLYRP4a SNV status, and the correlation between Faecalibacterium prausnitzi and Dopaquinone and 2ʹ-Deoxyguanosine.(C) The abundance of Bacteroides ovatus related to SNV status of LILRA6 and the correlation between Bacteroides ovatus and Desaminotyrosine, Tyrosyl-Leucine, and 3-(Uracil-1- yl)-L- alanine.

Three metabolites in the putative causal relationship network were involved in tryptophan and tyrosine metabolism, i.e. Tyrosyl-Leucine, Dopaquinone, and Desaminotyrosine. The level of Dopaquinone was causally related to the variation of PGLYRP4 (rs148195147, c.G602A, p.R201Q) and the abundance of Faecalibacterium prausnitzii (Figure 5(b)). Also, the Faecalibacterium prausnitzii is associated with the level of the 2ʹ-Deoxyguanosine, a metabolite involved in the folate biosynthesis pathway. A variation located in LILRA6 (rs56257556, c.A931 G, p.N311D) was associated with Desaminotyrosine, Tyrosyl-Leucine, and an Alanine metabolites 3-(Uracil-1-yl)-L-alanine mediated by Bacteroides ovatus (Figure 5(c)). LILRA6 is a leukocyte immunoglobulin-like receptor and plays roles in the pathways of the innate immune system and class I MHC mediated antigen processing and presentation. The de novo mutation of LILRA6 was previously observed in ASD.⁵⁸

Association of cytokines with host genetics and microbiota

Emerging research suggests that immune dysfunction is a risk factor contributing to the neurodevelopmental defects observed in ASD.⁵⁹ We then tested the levels of seven cytokines in the blood serum (Table S6). The level of IL-1β, IL-4, IL-17α, IL-10, and TGF-β was significantly higher in ASD children compared with TDs (Figure 6(a)), consistent with the previous studies.^8,17,18 Several studies have suggested that host genetic can reshape the gut microbiome and its metabolites through aberrant immune activation.⁸ The association between genetic variation and the differing cytokines was tested. The level of all the five differentially expressed cytokines was associated with SNVs (Figure 6(b), Table S7). Including known ASD-associated genes MUC2, MUC4, SETD5, MUC5B, INO80D, SEPT9, PGLYRP4, DST, and CHD6, among others. Besides, all the five cytokines were correlated with microbes and metabolites (Figure 6(c-d), Table S8-9). These results indicated that the interaction of gut microbes-metabolites and host genetics might play essential roles in the aberrant immune activation in ASD children.

Figure 6.

Cytokines in ASD and their association with gene, microbes, and metabolites

(A) The level of in cytokines in ASD and control samples, Wilcoxon test, ***fdr<0.001, **fdr<0.01, and *fdr<0.05.(B) The heatmap of associations between SNV (gene) and cytokines, * linear regression, p-value<0.05, ** p-value<0.01.(C) The heatmap of associations between microbes and cytokines, Spearman correlation, * p-value<0.05, ** p-value<0.01.(d) The heatmap of associations between metabolites and cytokines, Spearman correlation, * p-value<0.05, ** p-value<0.01.

Ethics statement

This study was approved by the Ethics Committee of Affiliated Yixing Hospital of Jiangsu University (Ethics NO, 2016055). All participants and their legal guardians were provided a written informed consent upon enrollment. Once the consent forms were signed, we screened them for eligibility criteria and sent questionnaires and sample collection kits to participants who meet the study eligibility criteria.

Study subject recruitment

The children with ASD were recruited as we described in our previous study⁷³. Briefly, ASD children were diagnosed according to the Diagnostic and Statistical Manual of Mental Disorders, 5th Edition. Each participant underwent physical, neurological, and behavioral examinations. Children were excluded from the study if they were previously diagnosed with a genetic condition such as tuberous sclerosis, significant sensory impairment, clinically evident inflammatory conditions, celiac disease, and other physiological conditions (i.e. depressive disorder, schizophrenia, and bipolar disorder). Participants that are on anti-inflammatory or anti-oxidant drugs, or take antibiotics, probiotics, prebiotics, and antifungal medications in 3 months before the feces collection were excluded. The Rome IV criteria for functional constipation were used for evaluating GI problems. The typically developed children (TD) for the fecal metabolism and cytokines measure were recruited from kindergartens through pediatricians and under physical, neurological, and behavioral examinations as the ASD group, subjects with any pathological conditions, including GI problems were excluded (Table S10).

Exon sequencing and data process

Peripheral blood samples were collected from the ASD and TD group for extracting DNA. Exome capture was performed using the Agilent SureSelect Human All Exon V6 r2 kit followed by Illumina paired-end sequencing.

The raw reads were quality-controlled using FastQC v0.11.2 and then mapped to the human reference genome (hg19) using bwa v0.7.17.⁷⁴ PCR duplicates were removed using Picard tools v2.18.12. SNVs were identified using FreeBayes v1.2.0⁷⁵ and filtered using bcftools v1.9. SNVs with base quality < 30, coverage < 30, with less than 3 reads mapped to the alternative alleles or located within 5bp nearby a gap were removed.

SNVs presented in the TD groups were removed. Then, SNVs were annotated by the annovar software to extract allele frequencies from the public database Allele frequency of these SNVs in the normal population was extracted from the 1000 Genome project,⁷⁶ the Exome Aggregation Consortium (ExAC) database,⁷⁷ and the Chines MillionomeI database (CMDB),⁷⁸ a large-scale Chinese genomics database including comprehensive variation and their allele frequency information from 141,431 unrelated healthy Chinese individuals. SNVs with MAF < 5% in these public databases and detected in more than two ASD samples were considered as candidates. GO enrichment analysis of genes with SNVs was conducted in R package.

Fecal sample collection

Fecal samples of each participant were obtained at the hospital or home and collected as directed in the instruction provided and frozen immediately until shipment. Samples were shipped using dry ice overnight to Nanjing Medical University and stored at −80°C until extraction.

Metagenomic sequencing and preprocess

About 2 μg DNA per sample was prepared for library construction and then sequenced on the Illumina HiseqX platform. Low quality or contaminated human reads were removed, and the high-quality reads were assembled into contigs using SOAPdenovo (v2.04).⁷⁹ MetaGeneMark (v2.10)⁸⁰ was used to predict genes from the assembled contigs. CD-HIT v4.5.8⁸¹ was used to generate a non-redundant gene catalog. Abundances of the genes were computed by aligning high-quality reads to the reference gene database. For taxonomic identity and functional assignment of unigenes, reads were aligned to the NCBI NR database (e-value≤1e-5) using DIAMOND (v0.9.9).⁸² The LCA algorithm⁸³ was used to conduct annotation. The detailed methods and results have been presented by our other study.

Gene functional annotations

The unigenes were functionally annotated by mapping to different functional protein databases with BLAST software. Predicted unigenes were assigned to the Kyoto encyclopedia of genes and genomes (KEGG), the evolutionary genealogy of genes: non-supervised orthologous groups (eggNOG), and carbohydrate-active enzymes database (CAZy) database using DIAMOND (v0.9.9). The abundances of each functional annotation were the sum of the abundance of annotation of each functional level.

Metabonomics analysis and data process

For each ASD and TDs, about 50 mg feces were mixed with 1000 μL of extract solvent, and samples were transferred to UHPLC-QE Orbitrap/MS analysis. The differential analysis was performed using the ropls package(v1.12.0), and metabolites were with VIP > 1, P-value (T-test) < 0.05 and FC>1.5 were defined as significant different in abundance. The detailed methods have been presented in our previous study.

Association analysis

Association test was performed using standard linear regression after adjusting for established confounding factors (age and sex). p-values were then corrected for multiple tests using the Benjamini-Hochberg method.

Causal inference analysis

The SNV-microbe-metabolite trios were assessed using the causal inference test (CIT)⁵⁷ to test the possibility that SNV causally influences microbe abundance and then the microbe metabolite. The enriched SNVs and their associated microbes, and the differential metabolites between ASD and TD group were subjected to the test. Briefly, the CIT has statistical tests for four conditions, all of which must be met for the microbe -mediated causal classification: (1) SNV and microbe are associated; (2) SNV is associated with microbe after adjusting for metabolite; (3) microbe is associated with metabolite after adjusting for SNV; and (4) SNV is independent of metabolite after adjusting for microbe. The CIT p-value was defined as the maximum of the component test p-values, and multivariate linear regression was used in the four-component tests.

Cytokine assay in healthy control and ASD patients

We measured cytokines in blood serum of the ASD cohort and an independent TD group including IL-1ß, IL-4, IL-6, IL-10, IL-17α, MCP-1, and TCG-β using BD™ CBA Flex Set (Table S6). Then the level of them was log-transformed, and the Mann-Whitney U test was applied to test the difference between ASD and TDs. These p-values were then corrected for multiple tests using the Benjamini-Hochberg method and cytokines with FDR<0.05 were considered as differentially expressed.