We recruited individuals with ASD and their parents, and unrelated neurotypical individuals as previously described.^35,36 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.

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.^37–39 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.

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).

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).

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: ASD=scale[log2(feature intensity+1)]+age+sex+batch 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.

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^45–47 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^49,50 and mzCloud (mzcloud.org)—to provide level 2 confidence (probable structure).⁴⁸

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: log2[(feature intensity+1)]=age+sex+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 (r_s) 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.

We estimated five r_s 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).⁵¹

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).

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 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).

We show the distribution of the r_ss within households in Figure 2. For HILIC data, mean r_s 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 r_s greater than 0.5. Proband-mother had a uniform-like distribution of r_s, 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 r_s 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 r_ss within households was similar to HILIC. Only five of all the features had r_s 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.

We show the within household r_ss 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 r_ss 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).

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.

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 r_ss (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 r_ss 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 r_ss were 0.31 and 0.68 for serum perfluorinated compounds and polybrominated diphenyl ethers, respectively.^57,58 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 r_s 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.²⁸,^59–62 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.⁶⁰