Plasma metabolomics of autism spectrum disorder and influence of shared components in proband families.

Prevalence of autism spectrum disorder (ASD) has been increasing in the United States in the past decades. The exact mechanisms remain enigmatic, and diagnosis of the disease still relies primarily on assessment of behavior. We first used a case-control design (75 idiopathic cases and 29 controls, enrolled at Boston Children's Hospital from 2007-2012) to identify plasma biomarkers of ASD through a metabolome-wide association study approach. Then we leveraged a family-based design (31 families) to investigate the influence of shared genetic and environmental components on the autism-associated features. Using untargeted high-resolution mass spectrometry metabolomics platforms, we detected 19 184 features. Of these, 191 were associated with ASD (false discovery rate < 0.05). We putatively annotated 30 features that had an odds ratio (OR) between <0.01 and 5.84. An identified endogenous metabolite, O-phosphotyrosine, was associated with an extremely low autism odds (OR 0.17; 95% confidence interval 0.06-0.39). We also found that glutathione metabolism was associated with ASD (P = 0.048). Correlations of the significant features between proband and parents were low (median = 0.09). Of the 30 annotated features, the median correlations within families (proband-parents) were -0.15 and 0.24 for the endogenous and exogenous metabolites, respectively. We hypothesize that, without feature identification, family-based correlation analysis of autism-associated features can be an alternative way to assist the prioritization of potentially diagnostic features. A panel of ASD diagnostic metabolic markers with high specificity could be derived upon further studies.

Introduction

Autism spectrum disorder (ASD) is a group of neurodevelopmental disorders characterized by impairment in communication and social interaction, and repetitive behaviors. Over the past decades, the prevalence of ASD in the United States has increased to 1 in 42 boys and 1 in 189 girls , with earlier, and inclusive, diagnosis contributing 12% and 56%–60% to the rise in prevalence, respectively. Since the relative stability of genetic factors cannot explain such a dramatic change in prevalence over such a short period, the steady increase may be partially due to environmental factors or interplay between gene and environment. Recent twin studies have suggested that nongenetic contribution of ASD could be between 30% and 50%. At present, a diagnostic test for ASD is not available through any objective biofluid test but rely on assessing children's behavior and developmental delays by physicians and other health specialists.,

The metabolome includes all of the small molecules found within a biological sample. The majority of the molecules are endogenous metabolites that participate in the essential chemical reactions, while exogenous chemicals that originate from the external environment can also be identified. Untargeted high-resolution mass spectrometry (HRMS) provides broad-spectrum coverage of the metabolome and therefore is commonly used to profile the molecular signatures presenting physiological status of an individual. Data-driven studies of the plasma metabolome have been successfully used to answer biomedical questions, for example, characterizing the disrupted physiological processes and molecular fingerprint of disorders and/or discovering the exogenous exposures associated with health outcomes.

Untargeted profiling of the metabolome has been applied to study ASD. These studies mostly employed a case–control design and aimed to identify biomarkers for early diagnosis, treatment, or prevent, and understand the pathophysiologic networks of ASD. To measure the blood or urine metabolome, analytical platforms such as gas chromatography–mass spectrometry, liquid chromatography–mass spectrometry, capillary electrophoresis–mass spectrometry, and nuclear magnetic resonance spectroscopy can be used. Each of these platforms has its own strengths and limitations in terms of speed, cost, sensitivity, and interpretability of the signals., and the final choice is dependent on the resources, sample types, and questions of interest in the studies. For exploratory analysis, employing multiple mass spectrometry platforms can maximize sensitivity and hence coverage of the metabolome, however, more than half of the published ASD metabolomics studies were using only a single platform.

Although powerful, untargeted metabolic profiling with HRMS has a major drawback—only a small fraction of the detected signals that have higher abundance in biofluids, such as glycolic acid, glycine, and arginine, can be routinely identified. Recent developments in mass spectral molecular networking and community-based spectral library have enhanced the interpretability of the metabolic profiles by providing correlation-based putative annotation, however, many of the features could still remain unknown. Furthermore, data-driven ASD studies typically generate a list of ASD-associated features and these potentially diagnostic features require structural confirmation and validation with independent cohorts. Since a majority of the pilot and exploratory studies employed a case–control design, associated features could be a mix of endogenous and exogenous molecules. Clearly, environmental chemicals have little biological relevance for reliable ASD prognosis or diagnosis, and the ability to distinguish the origins of features prior to tedious follow-up analyses would improve feature prioritization and accelerate pace of discovery.

We set two aims for this study. In Aim 1, we sought to discover new insights towards the molecular mechanism behind children with ASD using orthogonal liquid chromatography–mass spectrometry method to maximize the detection of plasma metabolome. In Aim 2, we posited that by studying the influence of shared genetic and environmental components in ASD families, we could deduce the origin and the role of the unknown ASD characterizing features. Specifically, we applied a metabolome-wide association study (MWAS) to an ASD-control comparison and conducted a series of family-based correlation analyses within and between probands and their parents. The overarching goal is to increase our understanding toward the etiology, physiology, and putative molecular biomarkers of ASD.

Methods

Study population

We recruited individuals with ASD and their parents, and unrelated neurotypical individuals as previously described., The MWAS used a case–control study design. We enrolled a total of 104 individuals (75 cases and 29 controls), who were part of a blood transcriptome study conducted at the Boston Children’s Hospital (BCH) in the years 2007-2012. We recruited unrelated controls from well-visit children to the Primary Care Center of BCH and from healthy individuals who visited the Division of Endocrinology of BCH for the evaluation of short stature. In our family-based correlation investigation, we enrolled parents of probands—40 mothers and 39 fathers. There were 31 families with at least one proband and at least one parent. BCH Institutional Review Board approved the study and we obtained informed consent from all participating subjects or their custodians prior to any data and sample collection.

Sample preparation and measurement

The whole blood was collected from participants in ethylenediaminetetraacetic acid-treated tubes and centrifuged at 2000 ×g for 10 min at room temperature to obtain plasma. These samples were then stored at −80°C until aliquoting and shipping. Plasma samples were thawed overnight in 4 °C refrigerator and then aliquoted into 0.5 mL screw-cap tubes prior to shipping in a dry iced box. Since blood samples were collected as part of standard clinical practice, time since last caloric intake varied from 5 min to 8 h. We prepared and analyzed the samples using established methods. We used dual-column chromatography (hydrophilic interaction liquid chromatography [HILIC; XBridge BEH Amide XP HILIC column; Waters, Waltham, MA, 50 × 2.1 mm, 2.5 μm] with positive electrospray ionization [ESI] and reversed-phase liquid chromatography [RPLC; C18 column; Higgins Analytical, Mountain View, CA, 50 × 2.1 mm, 2.6 μm] with negative ESI) coupled with HRMS for measuring the metabolome.

Data preprocessing

We converted raw data to open mzxml format using Proteowizard and conducted data preprocessing including peak detection, noise filtering, peak quantification and alignment, averaging signals of triplicates, peak matching and batch effect correction with apLCMS and xMSanalyzer. The pipeline generated a data table for each of HILIC and RPLC that comprises the relative signal intensities of 13 098 and 10 522 features, respectively.

We excluded features in the highest 20% of the coefficient of variations across quality control samples and those with greater than 80% zero intensity across study samples. The number of features available in our analyses was 10 173 (HILIC) and 8011 (RPLC). We presented the details about inclusion and exclusion criteria for ASD cases, controls and proband families, metabolome measurement, and analyses of quality control samples in the Supplementary material (Methods S1-S3).

Statistical analyses

Analytical overview

We divided the analysis into two parts as shown in Figure 1. (A) For Aim 1: metabolomic profiling of ASD and (B) for Aim 2: feature correlations in the shared environment. We conducted the analyses separately for HILIC and RPLC data. Using a case–control design, (A1) we executed an MWAS to unbiasedly discover features associated with ASD. Then, we conducted a series of analyses based on the MWAS selected features. (A2) We identified their chemical names through matching the measured mass-to-charge ratio (m/z) with the theoretical values. We leveraged a family-based design to understand the influence of shared genetic and environmental components on ASD. Specifically, we (B1) computed the correlation of the MWAS features in proband-mother and proband-father pairs and (B2) visualized the correlation patterns as correlation globes. Additionally, we provided complementary pathway enrichment for Aim 1 in the Supplementary material (Method S4).

Figure 1.

Overview of the data collection and statistical analyses conducted in this study. After collecting plasma samples of ASD cases, controls, and parents of ASD, we used two analytical platforms, HILIC and RPLC that were coupled with high-resolution mass spectrometry (HRMS) to conduct untargeted metabolic profiling. Raw data was processed, including feature detection, peak alignment, batch effect correction, and feature filtering before the data sets were sent to downstream analysis. For physiological analysis, we first run (A1) an MWAS on cases and controls to identify significant features associated with ASD (FDR = 0.05). Then we (A2) putatively named the compounds by matching their accurate masses to those unambiguously found in metabolite databases. Using the MWAS selected features, we conducted the shared environment analyses after adjusting for confounders (i.e., analyzed the extracted residuals from the adjusting model). We (B1) first investigated the Spearman’s rank correlations within households (i.e., in proband-mother, proband-father). Then we (B2) visualized the correlation patterns within households as correlation globes.

Metabolome-wide association study

To conduct an MWAS, we first extracted 75 cases and 29 controls from the data table. We filtered out features with zero variance in each group and analyzed the features that had greater than zero variance in both cases and controls. The number of available features was 10 143 (HILIC; 99.8% of the input features) and 7997 (RPLC; 99.8% of the input features). Then we scaled the log-transformed features and ran a generalized linear model as follows: where ASD is a binary variable denoting the diagnosis of ASD; age is a continuous variable (months); sex is a binary variable; batch is a categorical variable denoting the analytical batch identifier of the samples. We used false discovery rate (FDR) to correct for multiplicity when determining statistical significance.

Annotation of features

We used xMSannotator to annotate detected features. For each feature, the accurate mass was searched against those chemicals listed in Human Metabolome Database (HMDB), Kyoto Encyclopedia of Genes and Genomes (KEGG), and LIPID MAPS at a 5 ppm error-tolerant window. To avoid generating excessive false annotations, we only used hydrogen adduct for m/z matching (HILIC: M + H; RPLC: M-H). A single feature could be unambiguously matched with unique chemical identity or matched with multiple chemicals in the same or across different databases. We only reported the MWAS significant features with unambiguous putative annotation in the main table and provided the full putative annotation results in the Supplementary material (Table S1: HILIC; Table S2: RPLC). We followed the identification confidence levels in high resolution mass spectrometric analysis by Schymanski et al. (level 1: confirmed structure, level 2: probable structure, level 3: tentative candidate(s), level 4: unequivocal molecular formula, and level 5: exact mass) and reported all the metabolite identification with a level 5 confidence level unless otherwise specified. For the ASD-associated features that are believed to be endogenous metabolites based on the putative annotation, we further validated their identities as follows: we identified spectral peaks in the in-house pooled plasma samples with sufficient intensity for tandem mass spectrometry. Tandem mass spectrometry (MS/MS) fragmentation was then performed on the matched parent m/z and the resulting fragmentation MS/MS spectra were compared to those found in two public databases—MetFrag, and mzCloud (mzcloud.org)—to provide level 2 confidence (probable structure).

Correlations within households

We studied the household correlations separately for mother and father instead of averaging their signals to a single value (i.e., proband-mother and proband-father pairs). In each of the proband, mother, and father groups, we first filtered MWAS features that had > 80% zero intensity and fit a model to subtract linear contributions of age, sex, and batch:

After matching with family ID, we had 35 pairs from 31 families and 79 mutual features in each group (HILIC; 45 mutual features for RPLC data). The model residuals were used as the input for all the familial correlation analyses. We estimated the Spearman’s rank correlations (rs) of the mutual features in proband-mother and proband-father and corrected for multiplicity with FDR. We used the proportion of significant correlations as a metric to gauge the influence of shared components.

Correlation patterns within household

We estimated five rs matrices—proband, mother, father, proband-mother, and proband-father—using the residuals estimated from the previous step. In addition, using Euclidean distance, we performed hierarchical clustering analysis on the features of ASD and arbitrarily set to create 9 feature groups to aid pattern inspection. The group membership was applied to mother and father and features were sorted in the same order to aid visual inspection of the correlation patterns within households. We created the correlation globes with the R package circlize (v 0.4.5).

Computational settings

Our analyses took the following default settings unless otherwise specified. To control spurious findings due to multiple testing, we used a 5% FDR and reported the Q value, which is an FDR adjusted P-value. FDR was calculated separately for HILIC and RPLC analyses. We transformed, or scaled, each variable feaure to have mean zero and unit variance. We used a fudge factor of 1 (i.e., x + 1) to avoid directly log-transforming variables with zero intensity. Correlation refers to the estimation of Spearman’s rank correlations. We executed all analyses using the computing environment R (v 3.5.1). For reproducibility purposes, all analytic code is publicly available on GitHub via an MIT license (github.com/jakemkc/autism_shared_env).

Results

Demographic summary of the study population is shown in Table 1. The majority of the enrolled subjects were Caucasian and non-Hispanic with less than 5% Asian. We enrolled more individuals with ASD (75 versus 29) that were younger (98 versus 147 months) and with more male (83% versus 62%) when compared with neurotypical controls. Demographic characteristics of mother and father, including sample size (40 versus 39), age in months (493 versus 502), and ethnicity (98% versus 97% white) were similar.

Table 1.

Demographic characteristics of autism cases, parents of cases, and controls in this study

Characteristic	Control	Autism	Mother	Father
Number	29	75	40	39
Age, median (interquartile range), months	152 (34)	83 (51)	501.5 (94.7)	504 (63)
Male, number (%)	18 (62)	62 (83)	0 (0)	39 (100)
Ethnicity, number (%)
White	27 (93)	71 (95)	39 (98)	38 (97)
Asian	1 (3)	2 (3)	0 (0)	0 (0)
Unknown	1 (3)	2 (3)	1 (3)	1 (3)

Molecular profile of autism spectrum disorder

Table 2 shows the putatively annotated chemicals together with their test statistics. Out of the MWAS significant features (HILIC: 125, RPLC: 66), we could annotate 20 (16%) and 10 (15.2%) of them by their accurate mass for HILIC and RPLC, respectively. Only a few of these were associated with an increased odd of ASD (HILIC: 5/20; RPLC: 0/10). ORs for HILIC data spanned from less than 0.01 to 5.84 whereas those for RPLC were from less than 0.01 to 0.29. Annotated chemicals came from diverse chemical classes, such as piperidines, flavonoids, carboxylic acids and derivatives for HILIC data, and carboxylic acids, naphthofurans, and glycerolipids for RPLC data. Using the MWAS significant features to conduct pathway enrichment analysis, we found that only glutathione metabolism was affected (Table S4; P-value = 0.048).

Table 2.

Putatively identified metabolites from the MWAS of autism cases and controls and the correlations of metabolites within the ASD families using two different analytical platforms

Platform	Monoisotopic mass	OR (95% CI)a		Q value	dbIDb	Name	Source	Spearman correlationc
HILIC								ASD-Mother	ASD-Father
	129.1518	< 0.01	(< 0.01-0.03)	0.04	C01740	Octylamine; N-Octylamine; Monoctylamine	Exogenous	−0.186	0.276
	155.1310	<0.01	(< 0.01-0.01)	0.01	C06184	N-Methylpelletierine	Exogenous	0.364	0.255
	149.9987	0.01	(< 0.01-0.07)	0.04	HMDB42032	Thiodiacetic acid	Endogenous	−0.282	0.178
	468.9511	0.04	(0.01-0.15)	0.01	HMDB60640	Lamivudine-triphosphate	Endogenous	NA	NA
	85.0891	0.06	(0.01-0.19)	0.01	HMDB34301	Piperidine	Exogenous	−0.145	−0.142
	323.3188	0.10	(0.03-0.22)	0.01	HMDB34373	N-(14-Methylhexadecanoyl)pyrrolidine	Exogenous	NA	NA
	213.9879	0.10	(0.03-0.25)	0.01	HMDB06801	2-Oxo-3-hydroxy-4-phosphobutanoic acid	Endogenous	NA	NA
	141.0112	0.12	(0.04-0.28)	0.01	HMDB60688	Nornitrogen mustard	Endogenous	NA	NA
	193.0600	0.13	(0.04-0.31)	0.02	C16789	Toxoflavine	Exogenous	NA	NA
	350.0096	0.15	(0.05-0.34)	0.02	HMDB37851	Apigenin 7-sulfate	Exogenous	0.346	0.299
	261.0402	0.17	(0.06-0.39)	0.02	HMDB06049	O-Phosphotyrosine	Endogenous	NA	NA
	341.9739	0.19	(0.08-0.39)	0.01	C18968	Carbophenothion	Exogenous	0.308	0.396
	276.0124	0.24	(0.11-0.46)	0.02	C11570	2-(2-Chloro-phenyl)-5-(5-methylthiophen-2-yl)-134oxadiazole	NA	NA	NA
	204.9810	0.26	(0.13-0.49)	0.02	C12284	Saccharin sodium anhydrous	Exogenous	NA	NA
	378.1943	0.30	(0.16-0.54)	0.02	HMDB61127	4R-Hydroxy solifenacin	Endogenous	NA	NA
	355.9497	0.31	(0.15-0.58)	0.04	HMDB15610	Silver sulfadiazine	Exogenous	0.204	0.331
	263.0844	2.89	(1.66-5.73)	0.05	HMDB42056	Tulobuterol	Exogenous	0.023	0.236
	500.3866	3.14	(1.74-6.37)	0.04	C15379	3beta-Hydroxylanostane-711-dione acetate	Exogenous	−0.020	0.078
	471.0923	3.49	(1.82-7.79)	0.05	HMDB61083	desbutyl-lumefantrine	Endogenous	−0.185	−0.116
	650.2211	5.84	(2.50-17.5)	0.03	LMPK12050358	5-Hydroxy-734-trimethoxy-8-methylisoflavone 5-O-neohesperidoside	Exogenous	0.241	0.195
RPLP
	326.2093	<0.01	(<0.01-<0.01)	0.05	HMDB31963	(3b6b8a12a)-812-Epoxy-7(11)-eremophilene-6812-trimethoxy-3-ol	Exogenous	0.003	0.312
	230.0579	<0.01	(<0.01-<0.01)	0.02	LMPK13110003	Visnagin	Exogenous	−0.115	−0.037
	410.1818	0.06	(0.01-0.18)	0.02	C08093	Oseltamivir phosphate	Exogenous	NA	NA
	274.2144	0.11	(0.03-0.28)	0.02	C13854	1-Dodecanoyl-sn-glycerol	NA	−0.155	−0.074
	217.0773	0.17	(0.06-0.36)	0.02	HMDB15328	Captopril	Exogenous	0.211	0.466
	216.0746	0.19	(0.08-0.38)	0.02	C17359	8-Hydroxyalanylclavam	Exogenous	0.281	0.498
	269.1991	0.21	(0.09-0.40)	0.02	HMDB32255	N-(Ethoxycarbonyl)methyl)-p-menthane-3-carboxamide	Exogenous	NA	NA
	322.1780	0.24	(0.11-0.48)	0.03	HMDB32702	Zeranol	Exogenous	NA	NA
	386.2305	0.25	(0.11-0.48)	0.03	C11990	Oleandolide	Exogenous	NA	NA
	332.2563	0.29	(0.15-0.53)	0.03	C19621	Floionolic acid	Exogenous	NA	NA

Only features with level 5 identification confidence (exact mass) and unambiguous compound matching with Human Metabolome Database (HMDB), Kyoto Encyclopedia of Genes and Genomes (KEGG), and LIPID MAPS are shown. For silver sulfadiazine, it is now located in the Toxin and Toxin Target Database (T3DB) with an accession number T3D3068. All of the metabolites were matched by exact mass through comparison with records in reference databases. Only O-phosphotyrosine was identified to level 2 (probable structure) by matching the MS/MS spectrum. A putatively annotated metabolite by exact mass, 2-oxophytanic acid, was not shown because its tandem mass spectrometry (MS/MS) spectrum did not match with the records in databases. Metabolites are sorted by size of the OR.

Input features to the model were log-transformed and scaled to have mean zero and unit variance.

Unique ID for the corresponding chemicals found in HMD, KEGG or LIPID MAPS.

Spearman’s rank correlations of annotated metabolites within the ASD families (i.e., in proband-mother, proband-father). NA represents the metabolite not consistently detected within the families.

Correlations in the shared environment

We show the distribution of the rss within households in Figure 2. For HILIC data, mean rs was 0.08 (range: −0.35 to 0.71) and 0.11 (range: −0.28 to 0.60) for proband-mother and proband-father, respectively. Only a few of the features had rs greater than 0.5. Proband-mother had a uniform-like distribution of rs, whereas that for proband-father is more resembled to normal. Of the 79 investigated features, 5% were significantly correlated—two were significant in proband-mother and proband-father, respectively, after FDR correction (i.e., four non-overlapping features). For RPLC data, mean rs was 0.12 (range: −0.23 to 0.80) and 0.13 (−0.29 to 0.58) for proband-mother and proband-father, respectively. Distribution of rss within households was similar to HILIC. Only five of all the features had rs greater than 0.5. For all 45 features investigated, 13% were significantly correlated—one was significant in proband-mother and five were significant in proband-father after FDR correction. All the significant correlations were positively correlated with a magnitude greater than 0.45.

Figure 2.

Violin plots showing the distributions of the correlations of ASD–associating features within the ASD families. (A) HILIC platform; (B) RPLC platform. Using the significant features found in the MWAS, we estimated the Spearman’s rank correlation for each feature between proband-mother, proband-father. The number of shared features for estimating correlations were 79 (HILIC) and 45 (RPLC) after data preprocessing for familial analysis. In each plot, we overlaid with a one-dimensional scatter plot and a boxplot showing median, interquartile ranges, and whiskers extending to the largest values within 1.5*interquartile range.

We show the within household rss of the MWAS significant metabolites in Table 2. Overall, only two out of seven (28.6%) endogenous and 13 out of 21 (61.9%) exogenous metabolites were consistently detected in the ASD families. The median rss within the families (across HILIC, RPLC, ASD-mother, and ASD-father) were −0.15 and 0.24 for endogenous and exogenous metabolites, respectively.

We show the within household correlation patterns in Figure 3. In the correlation globes depicting the patterns of 79 features from HILIC data, we did not observe strong similarity in the correlation patterns between proband and mother (Figure 3A) nor between proband and father (Figure 3B). Instead, we observed general similarity in patterns between mother (Figure 3A, right half) and father (Figure 3B, right half), and to a lesser extent, between proband-mother (Figure 3A, across the globe) and proband-father (Figure 3B, across the globe). These observations of correlation patterns were also found for the 45 features from RPLC data (Figure 3C and 3D).

Figure 3.

Correlation globe showing the correlation patterns in the shared environment. (A) and (B): HILIC platform; (C) and (D) RPLC platform. Using the HILIC platform as an example, we estimated the Spearman’s rank correlations of 79 ASD correlating features shared in ASD families (RPLC: 45 shared features). To facilitate visual inspection of the patterns, we assigned the features in ASD, mother, and father into nine feature groups based on the results from hierarchical clustering on ASD cases. Each correlation globe is showing the correlations within ASD, within parents (mother or father), and between ASD and parents (mother or father). Features are arranged as a circular track. Left and right halves of the globe represent features in ASD and parent (mother or father), respectively. Only Spearman’s rank correlations greater than 0.5 and smaller than −0.5 are shown as connections in the globe. Red line denotes positive correlation, and dark green line denotes a negative one. Color intensity and line width are proportional to the size of the correlation. Within-group and between-group correlations are shown outside and inside of the track, respectively. Correlations between ASD and parent are indicated by the lines linking across the vertical half of the globe.

Discussion

Molecular profiling of ASD

In this study, we employed orthogonal separation techniques coupled with HRMS to maximize coverage of the human plasma metabolome. Our comprehensive analyses have shown that a total of 191 chemical features were associated with ASD in our study cohort. These markers could be endogenous molecules acting as mediators in the causal pathways for ASD, noncausal indicators of ASD, or exogenous exposures with natural or anthropogenic origins. In the pathway enrichment analysis, we found that glutathione metabolism was perturbed in ASD.

In our initial annotation, we found three putative MWAS significant metabolites detected using HILIC could be endogenous compounds (Table 2). Upon further identification, we found that 2-oxophytanic acid was a false positive and concentration of 2-oxo-3-hydroxy-4-phosphobutanoic acid was too low for MS/MS comparison. Only O-phosphotyrosine (pTyr) was identified by a matched mass spectrum in the reference database.

pTyr (OR 0.17, 95% confidence interval [CI] 0.06-0.39) is a product of tyrosine phosphorylation, which is the addition of a phosphate group to tyrosine and is catalyzed by tyrosine kinases. Phosphorylation of tyrosine residues in proteins is an important post-translational modification that is implicated in various biological processes including cell–cell signaling and proliferation as well as neuronal maturation and synaptic plasticity. Concentration of pTyr in body fluid could be correlated with tyrosine kinase and pTyr phosphatase activities in tissue. In the current study, sorting out the proteins with tyrosine residue responsible for lower concentration of pTyr in plasma samples from ASD warrants further studies. Nonetheless, a receptor tyrosine kinase encoded by MET is implicated in pathophysiological changes of ASD. A non-coding promoter variant of MET, rs1858830 C allele is reported as a genetic risk factor of ASD. rs1858830 CC genotype decreases MET promoter activity, which results in down-regulation of MET gene expression.

2-oxo-3-hydroxy-4-phosphobutanoic acid (OR 0.10, 95% CI 0.03-0.25) is part of the vitamin B6 pathway. Some have suggested that vitamin B6 can be used as a supplement in brain disorder. However, a randomized controlled study to investigate this relationship had shown negative results.

We found a handful of exogenous chemicals including pharmaceuticals, natural dietary molecules, and food additives to be significantly associated with ASD (Table S4). The combination of case–control design and blood sampling of ASD cases after diagnosis make our findings merely suggestive about their potential roles in the cause and development of ASD. Nonetheless, these shortlisted chemicals could be good candidates for future etiological research.

Familial correlations

The metabolome of ASD cases is affected by genetics, the environment, and the disease. We found 191 features were associated with ASD in our case–control MWAS analysis. Many of these features were not putatively annotated and we hypothesize that these could be either (1) exogenous chemicals from confounders such as unique environment or (2) endogenous biomarkers that are pathophysiological signatures of ASD and possess diagnostic value.

We leveraged a family-based study design and sought to obtain hints about the origins of the ASD-associated features. Key findings of our correlation analyses include low rss (averaged across HILIC and RPLC data: proband-mother; median [interquartile range], 0.09 [0.35]; proband-father; median [interquartile range], 0.08 [0.3]) and low similarity in correlation patterns in proband-mother and proband-father. In a familial setting, assuming family members stay 12 h a day at home, proband share about 50% of the genetics and the environment with the parents. In other studies, the median rss of exogenous chemicals between couples were about 0.21 (serum; polychlorinated biphenyls, organochlorine pesticides, polybrominated chemicals, per- and poly-fluoroalkyl substances, and metals) and 0.41 (urine; triclosan, phenols, bisphenols, parabens, and phthalates). For child–parent pairs, the median rss were 0.31 and 0.68 for serum perfluorinated compounds and polybrominated diphenyl ethers, respectively., In addition, almost all of the exogenous chemicals were positively correlated in these studies.

Further, using the known origin of some of the putatively identified metabolites (Table 2), we found that only 28.6% of the ASD-associated endogenous metabolites were also detected in parents, which is much smaller than that for exogenous metabolites (61.9%). Also, the median rs within the families for endogenous metabolites was smaller than that of exogenous metabolites (−0.15 and 0.24, respectively). These comparisons and analyses suggest that familial correlation in the ASD families could be used to infer the origin of the ASD-associated features found in the case–control MWAS. Correlation has a range between −1 and +1 and we hypothesize that the smaller the familial correlation of an MWAS feature, the more likely that it is endogenous, and may rank higher in priority for follow-ups such as structural confirmation and validation of diagnostic value.

Several studies have shown the potential of using plasma metabolites as clinical biomarkers for diagnosing neurological diseases., For example, West et al. found a set of over 100 metabolites that were able to discriminate ASD subjects from typically developing children with a maximum accuracy (area under the curve) of 81%. In an Alzheimer disease study, Stamate et al. reported that the diagnostic performance of using plasma metabolites had the potential to match the well-established cerebrospinal fluid biomarkers.

Limitations

We identified a few study limitations that pose a challenge to our findings. First, individuals with ASD were not matched with neurotypical controls by age and by sex; specifically, the control group was older with a higher proportion of female. To address this case–control imbalance, we adjusted the analyses by age and sex. Second, blood samples were not collected under an overnight fasting protocol. Our analysis may suffer from higher variability and hence reduced statistical power to detect associations. Third, for most of the putatively annotated metabolites, we did not conduct additional identification such as matching fragmentation patterns in a spectral library or comparison with authentic standards and thus they should be interpreted with caution. Fourth, autism is a complex and multifactorial disease and may require a larger sample size to capture a greater number of metabolites at smaller association sizes, that in totality, may ascribe more variation to autism. The sample size in this study (n = 104) is insufficient for a reliable subtype analysis. Last, our MWAS analysis could be confounded by unmeasured factors such as medication and dietary factors. The primary goal of this exploratory study is to investigate what deep metabolomics measurement could tell us about autism and the families with children who have autism.

Conclusions

Our case–control MWAS analysis revealed that a total of 191 features were associated with ASD. An identified metabolite, O-phosphotyrosine, was associated with a decreased risk of ASD. We also found glutathione metabolism was affected in ASD. Family-based correlation of ASD-associated features can assist the prioritization of potentially diagnostic features. Further studies are required to select a panel of reproducible metabolites, quantify their the diagnostic performance in different clinical settings of the metabolites, and identify the exogenous environmental factors associated with the metabolites.

Click here for additional data file.