Subject Samples

The experimental subjects were initially recruited through the UC Davis M.I.N.D. Institute Clinic, Regional Centers, referrals from clinicians, area school districts and community support groups such as Families for Early Autism Treatment (FEAT). Subjects were limited to an age range of 4–6 years. Typically developing participants were recruited from area school districts and community centers. All facets of this study were approved by the University of California at Davis Institutional Review Board (IRB). Written informed consent was obtained from the parent or guardian of each participant and data were analyzed without personal information identifiers then subjects completed diagnostic and psychological measures. Study participants with ASD were enrolled under inclusion criteria consisting of a diagnosis of autism spectrum disorder based on the DSM-IV criteria determined by an experienced neuropsychologist (BAC), which was further corroborated by the following measures using research reliable clinicians: the Autism Diagnostic Observation Schedule-Generic (ADOS-G) provides observation of a child's communication, reciprocal social interaction, and stereotyped behavior including an algorithm with cutoffs for autism and autism spectrum disorder; the Autism Diagnostic Interview-Revised (ADI-R) is a comprehensive, semi-structured parent interview that assesses a child's developmental history and relevant ASD characteristic behaviors and generates a diagnostic algorithm for children with ASD. Based on the DSM-IV criteria [21], only children with strictly defined autistic disorder were enrolled whereas children with pervasive developmental disorder-not otherwise specified (PDD-NOS) or Asperger Syndrome were excluded from the study. The Social Communication Questionnaire (SCQ) was used as a screening tool to ensure the absence of symptoms of ASD in the TD children. The patients recruited for this study were primarily Caucasian and the ages were similar between groups. However, the participants with autism had lower IQ scores than the TD subjects [22], [23].

The exclusion criteria for all subjects included the presence of Fragile X or other serious neurological (e.g., seizures), psychiatric (e.g., bipolar disorder) or known medical conditions such as autoimmune disease and inflammatory bowel diseases/celiac disease. All subjects were screened via parental interview for current and past physical illness. Children with known endocrine, cardiovascular, pulmonary, and liver or kidney disease were excluded from enrollment in the study. Dietary restriction for participation in the study was not required with the exception of an overnight fast. Participation in the study required two clinical visits for behavioral assessment and blood draws. After application of exclusion criteria, the final study group consisted of 104 children, 69 with ASD and 35 in the TD group.

Samples were collected on Thursday morning visits to the M.I.N.D. Institute over a period of 13 months. Blood was drawn into a 9.6 mL EDTA Vaccutainer tube by an experienced pediatric phlebotomist between the hours of 8 and 10 AM following an overnight fast. Tubes were immediately inverted 6–8 times to assure mixing with the anticoagulant and placed on ice. Immediately after plasma separation and aliquoting, samples were sent on the morning of the draw via courier with a barcode label, wrapped tube cap with a strip of parafilm; bubble wrapped then set in a biohazard bag which was placed inside a carrier between coolant packs. Samples were stored at −80°C.

Samples from 103 of the 104 children were sent to Stemina for metabolomic analyses on dry ice… Upon receipt, 5 samples were removed after visual inspection and observation of overt hemolysis and the remaining 98 samples were analyzed by mass spectrometry. Quality checks of the raw mass spectrometric data from the 98 samples were performed, resulting in removing data from 16 patient samples that did not contain MS data from all 5 methods from further analysis. The final 82 samples used in these studies originated from 52 children with ASD and 30 children in the TD group. The children were chosen so that the age and gender distributions were similar across the groups. There was no statistical difference in age between ASD cases and the TD children for the current study (Welch's t-test p = 0.25) (see Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Patient demographic information.

https://doi.org/10.1371/journal.pone.0112445.t001

Regarding patient medication, 18 out of 52 of the subjects with ASD in this study were taking medications which included Risperidone (5), Sertraline (3), Aripiprazole (2), antihistamines (2), antivirals (2), antifungals (2), and various other less frequent drugs. Three of the 30 typical subjects were taking medications, which included methylphenidate (1), albuterol (1) and loratadine (1). Ten of the 52 ASD subjects were on a gluten and/or casein-free (GFCF) diet. Importantly, blood draws were administered prior to eating and any morning administration of any medication.

Sample Preparation for LC-MS

Plasma samples were split into 50 µl aliquots and stored at −80°C prior to metabolite extraction. Samples were kept on ice during these procedures. Samples were randomized into three batches for the LC-HRMS analysis such that diagnosis, IQ, age and ethnicity were equally distributed in each batch. Small molecules were extracted from 50 µL plasma aliquots using 450 µL of 8∶1 methanol: water solution at −20°C [24]. The extraction solution also contained internal standards. The samples were agitated for 10 minutes at 2–8°C then centrifuged at 18,400×G for 20 minutes at 4°C to remove the precipitate. The supernatant was transferred to a fresh tube and the centrifugation step was repeated to remove any residual precipitate. After the final centrifugation, 450 µL of supernatant was transferred to a fresh tube then evaporated to dryness in a SpeedVac, then resolublized in 45 µL of a 50∶50 mixture of 0.1% formic acid in acetonitrile: 0.1% formic acid, also containing internal standards. This solution was then transferred to a high performance liquid chromatograph (HPLC) autosampler injection vial for LC-HRMS analysis.

Mass Spectrometry

Both targeted GC-MS as well as untargeted LC-HRMS were employed for better metabolome coverage. Four untargeted LC-HRMS methods were used including C8 or HILIC chromatography coupled to electrospray ionization in both positive and negative ion polarities, resulting in 4 separate data acquisitions per sample. For each methodology and condition, only a single sample aliquot was assessed, due to limited material availability.. LC-HRMS methods were developed and tested prior to the evaluation of the clinical patient samples to optimize the breadth of coverage of small molecule metabolites.

Liquid Chromatography High Resolution Mass Spectrometry

LC-HRMS was performed using an Agilent G6540 Quadrupole Time of Flight (QTOF) system consisting of an Agilent 1290 HPLC coupled to a high resolution (QTOF) mass spectrometer. Electrospray ionization (ESI) in both positive and negative ion modes was employed using a dual ESI source under high-resolution exact mass conditions. 2 µL of sample was injected. A Waters Acquity ultra high performance liquid chromatography (UPLC) BEH Amide column with dimensions 2.1×150 mm, 1.7 µM particle size was used for Hydrophilic Interaction Liquid Chromatography (HILIC), and maintained at 40°C. Data was acquired for each sample for 29 minutes at a flow rate of 0.5 mL/minute using a solvent gradient with 0.1% formic acid in water and 0.1% formic acid in acetonitrile. An Agilent Zorbax Eclipse Plus C8 2.1×100 mm, 1.8 µM particle size column was used for C8 chromatography and data was acquired for each sample for 50 minutes at a flow rate of 0.5 mL/minute using a gradient with 0.1% formic acid in water and 0.1% formic acid in acetonitrile and maintained at 40°C.

Gas Chromatography - Mass Spectrometry

GC-MS analyses were performed at the West Coast Metabolomics Center at UC Davis as described in [25]. GC-MS data was acquired using an Agilent 6890 gas chromatograph coupled to a LECO Pegasus IV TOF mass spectrometer. Metabolite identification was done by comparing sample data to a database of over 1,000 compounds identified by GC-MS that includes mass spectra, retention indices, structures and links to external metabolic databases.

Metabolite chemical structure confirmation by LC-HRMS-MS

The chemical structures of key metabolites were further confirmed using tandem mass spectrometry (LC-HRMS-MS) methods with chromatographic conditions identical to those used for their discovery. LC-HRMS-MS analyses were performed on an Agilent QTOF mass spectrometer for patient samples and/or, reference blood samples with collision energy conditions optimized to obtain the highest quality product ion spectra. The resulting product ion spectra were then compared to MS-MS spectra available in public spectral databases such as METLIN [26], MassBank [27] and Stemina's own SteminaMetDB database.

Data Analysis

Classification Modeling.

Model development was performed with two primary goals: 1) to robustly rank the importance of metabolites in discriminating ASD using a VIP (Variable Importance in the Projection) score index and 2) to identify the minimum set of predictive metabolites needed to reach the highest levels of differentiation of the ASD and TD experimental classes. The final models were created by training a Partial Least Squares Discriminant Analysis (PLS-DA) or Support Vector Machine (SVM) classifier using the entire 61-sample training set. The modeling techniques PLS-DA as well as SVM with a linear kernel [33], [34] were both utilized to demonstrate that the molecular signature can be predictive using multiple approaches. PLS and SVM classification models were created using the R package Classification and Regression Training “caret” version 5.17–7 [35]. Receiver operator Curve (ROC) analysis was performed using the R package ROCR version 1.0–5 [36].

A nested cross validation (CV) approach (Figure 1) was used to meet the first objective of model development - a robust measure of feature VIP scores. The 179 features from the 61 sample training set were analyzed using 100 resamples with an 80∶20 split to weight the importance of each of the 179 statistically significant features. The tuning loop utilized 10-fold cross validation to tune model parameters (cost parameter C for SVM and the number of components for PLS-DA). The recursive feature elimination loop was used to identify the best performing feature subset from each iteration using steps of 20 features. The results from the 100 resamples were used to estimate model performance and create a robust biomarker VIP score index to rank the importance of each of the 179 features in classification of ASD from TD individuals.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Classification modeling process.

A three-layer nested cross-validation approach was applied using both PLS-DA and SVM modeling methods to determine significant features capable of classifying children with ASD from TD children. The 179 features of the training set were analyzed using a leave-one-group-out cross-validation loop as described. The results from this cross-validation process were used to estimate model performance and create a robust feature VIP score index to rank the ASD vs TD classification importance of each of the 179 features. These feature ranks were used to evaluate the performance of the molecular signature using an independent validation set.

https://doi.org/10.1371/journal.pone.0112445.g001

Feature VIP robustness was measured by resampling the training set 100 times using an 80∶20 split into 49-sample CV training and 12-sample CV test sets. VIP scores were calculated for each of the 100 resamples and the most informative features at each resample were identified by backwards recursive feature elimination (in 20-feature steps) using the Area Under the ROC Curve (AUC). The most informative set of features was then used to predict each CV test set. The VIP scores were averaged across the 100 resamples to create the VIP index for each feature. The classification performance metrics of the CV test sets were averaged across resamples to understand potential future performance.

The second objective of the classification modeling approach was to identify the minimum number of features with the highest level of classification accuracy. This objective was met using feature subsets based on the ordered VIP score index and evaluating the subset performance in the validation set of samples. The classification models were created using the entire 61 sample training set and by stepping through features. The feature stepping process utilized the 20 top VIP features then added the next 20 highest weighted features until all 179 features were evaluated.

Performance metrics (Accuracy, Sensitivity, Specificity, and ROC analysis) were determined based on the prediction of the 21 sample independent validation set for assessment of the molecular signature at each feature subset bin size (Table S3). Accuracy is defined as the proportion of correctly classified participants and is calculated by dividing the number of correctly classified participants by the total number of participants in a sample set. Specificity is the proportion of correctly classified TD individuals out of all TD participants in a sample set. Sensitivity is the proportion of correctly classified ASD individuals out of all participants with ASD in a sample set. The top 179 features were also compared for rank between SVM and PLS modeling methods (Figure 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Feature Importance Rankings.

The top 179 features were compared for rank between SVM and PLS modeling methods. The lowest rank scores represent the most important features.

https://doi.org/10.1371/journal.pone.0112445.g002

Training Set Model Performance

SVM and PLS classification methods were used to discriminate between samples from children with ASD and TD children using the 179 selected features as variables and each feature's contribution toward classification was evaluated for future biomarker development efforts. Based on the best performing model from each of the 100 nested CV resampling iterations, ROC plots were generated for the average of the 100 resamples to understand performance of each modeling method (SVM and PLS-DA). Both SVM and PLS modeling methods indicated that a metabolic signature could be detected that could classify children with ASD from TD individuals (Table S2). For the 61-sample training set, the average ASD prediction accuracy of the SVM model was 0.86, with AUC values of 0.95 (95% CI 0.94–0.96). The PLS model gave an average prediction accuracy of 0.84 with AUC values of 0.92 (95% CI 0.91–0.94). To confirm that the model classification accuracies were not random results, the modeling process was repeated with random permutations of the diagnosis class labels. These results showed near random classification, with AUC values of 0.52 (95% CI 0.48–0.57) and 0.52 (95% CI 0.49–0.56) for SVM and PLS, respectively, indicating that the 179 features did not discriminate the classes by chance (Figure 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Performance of the SVM and PLS models.

Average AUC and accuracy of the (a) SVM and (b) PLS models containing different numbers of features. The bar graphs show the number of optimal models which were derived from recursive feature elimination process that was included in the resampling process for the indicated number of features.

https://doi.org/10.1371/journal.pone.0112445.g003

Anticipating that blood tests for ASD may be more efficient and less expensive if they measure an optimally lower number of metabolites, the classification modeling paradigm also included a feature number optimization in each model, based on the highest resulting AUC. The feature sets were evaluated using feature subsets based on the ordered VIP scores of individual features to identify the minimum number of features that maximized performance for each modeling method (Table S2). These data together indicate that not all of the features contributed equally to the models and that the number of features could be reduced by removing those that contributed less while still retaining model accuracy and robustness. As a result, the entire set of 179 features was not required for optimal model performance for either of the modeling methods (Figure 3). The results from the model training process indicated that SVM models that were trained using an 80 feature set exhibited the best combined classification performance metrics (when compared to PLS and other SVM results) with an average accuracy of 90%, an average sensitivity of 92%, an average specificity of 87%, and an average AUC of 0.95 (Table S2).

Validation Set Model Performance

Different subsets of features, created based on the weighted VIP scores, were evaluated independently of the outer cross-validation loop using the 21-sample independent validation set. The 80-feature SVM model described above had a classification prediction accuracy of 81%, a sensitivity of 85%, a specificity of 75% and an AUC of 84% (Figure 4, red line; Table 3). The best performing PLS model, comprised of 160 variables, had an accuracy of 81%, a sensitivity of 92%, a specificity of 63% and an AUC of 0.81% (Figure 4, blue line; Table 3). Detailed results are shown in Table S3. The results suggest that at least 40 features are needed to reach an accuracy of 70% and that a range of 80 to 160 features had the best performance with this independent validation sample set as well as the training set of samples.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. ROC curve performance of the classification models from the training and validation sets.

The average of 100 iterations of the classifier for the best performing feature sets following recursive feature elimination comparing ASD vs. TD samples (Black and Grey Lines). The blue (PLS) and red (SVM) lines are ROC curves of the best performing validation feature subsets. Vertical bars represent the standard error of the mean.

https://doi.org/10.1371/journal.pone.0112445.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Classifier performance metrics based on predictions on the independent 21-sample validation set, showing the feature sets with the highest accuracy.

https://doi.org/10.1371/journal.pone.0112445.t003

Confirmation of Metabolite Chemical Structures

The chemical identities of the 7 LC-MS mass features that were confirmed by LC-HRMS-MS are shown in Table 4. Included in the metabolites confirmed by LC-HRMS-MS or targeted GC-MS was homocitrulline, which was decreased in ASD patients. Other metabolites showing significant up or down regulation in our study include: aspartate, glutamate, DHEAS, citric acid, succinic acid, methylhexa-, tetra- and hepta-decanoic acids, isoleucine, glutaric acid, 3-aminoisobutyric acid and creatinine. These are listed in Table 4 and represent a variety of molecular classes including amino acids, organic acids, sterols, and fatty acids.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Confirmed metabolites.

https://doi.org/10.1371/journal.pone.0112445.t004

Identification of metabolites previously associated with ASD

Examples of metabolites showing significant up or down regulation in the current study that have been previously associated with autism include:

• Tricarboxylic acid cycle associated molecules including citric acid (decreased) and succinic acid (increased) were found to be significantly altered in the ASD participants. Elevations in urinary succinate [16], [17] and decreased urinary citrate [37] in children with autism have been previously reported.
• Fatty acids have previously been observed to be decreased in the plasma of children with ASD, similar to our observations for methylhexa-, tetra- and hepta-decanoic acids [12]. Links between saturated fatty acid metabolism and oxidative stress have been reported in erythrocytes in children with ASD [38].
• 3-aminoisobutyric acid was increased in samples from participants with ASD. This is also consistent with previous findings [39].
• Creatinine was decreased in children with ASD and is consistent with the findings of Whitely et al, who observed similar changes in urinary creatinine in children diagnosed with PDD [40].

Evidence for a role in mitochondrial dysfunction in ASD

The goal of this study was to evaluate biomarkers in blood. When metabolism is disrupted, active transporters and tissue specific differences in metabolism can cause different levels of the same metabolites in different biological compartments (plasma, urine, CSF, etc.). When discussing individual metabolites in the context of autism, it is important to recognize that autism is a systemic disease that affects other organ systems besides the brain. For example, serotonin levels have been reported to be elevated in blood in some patients with autism, and other evidence suggests that intracerebral serotonergic activity is decreased in ASD. Serotonin does not cross the blood brain barrier, and the brain has different enzymes for serotonin synthesis than peripheral tissues [41].

Many of the confirmed metabolites that are associated with ASD are relevant to aspects of mitochondrial biology. Mitochondrial disease or dysfunction may be a risk factor for autism [42]. In addition, several other observed metabolites are associated with other processes already proposed to be involved in ASD including oxidative stress [43] and energy production [44].

• Aspartate and glutamate levels in blood were significantly elevated in ASD, as has been observed in previous studies [45], [46]. Mutations in the aspartate/glutamate mitochondrial transporter, SLC25A12, have been previously associated with ASD [47]. This transporter is an important component of the malate/aspartate shuttle, a crucial system supporting oxidative phosphorylation, adenosine triphosphate production, and key metabolites for the urea cycle [47].
• DHEAS, the predominant plasma sterol, was increased in children with ASD. DHEA is known to affect mitochondrial energy production through inhibition of enzymes associated with the respiratory chain [48] with variable findings in children with ASD [49], [50].
• The branched chain amino acid isoleucine was reduced in samples from children with ASD versus TD children as observed by others [51]. Possible molecular mechanisms would include mutation in the branched chain amino acid kinase dehydrogenase (BCKD-kinase), a mitochondrial enzyme [52] as well as amino acids in energy metabolism [53].
• Glutaric acid levels were elevated in children with ASD. Increased urinary glutaric acid occurs in a variety of neuronal deficiencies such as glutaryl-CoA dehydrogenase (GCDH) deficiency. A significant portion of the glutaric acid metabolism takes place in the mitochondria [54].

The potential relationship of the gut microbiome with ASD

This potential connection between the gut microbiome and ASD is also receiving considerable attention [55]. Metabolomic studies of urine from individuals with ASD have identified molecules associated with the microbiome such as dimethylamine, hippuric acid and phenylacetylglutamine [16], [17]. We observed decreased plasma levels of p-hydroxyphenyllactate, a metabolite associated with bifidobacteria and lactobacilli that is known to serve as an antioxidant both in the circulation and tissues [56]. We have yet to identify other microbiome related metabolites.

Novel metabolic alterations in ASD

We identified novel statistically significant changes in some metabolites that had not been previously reported in other metabolomics profiling as well as identifying novel changes in metabolites that had never before been associated with ASD. Significant changes in the levels of aspartate, citrate, creatinine, DHEA-S, hydroxyphenyllactate, indoleacetate, isoleucine, glutamate and glutarate between ASD and TD individuals were identified in this study compared to previous metabolomics studies of urine metabolites where changes in these metabolites were not significant [17]. These differences could be related to transport and accumulation of metabolites in the urine compared to the blood, differences in the study populations, or application different LC-MS and GC-MS methodology allowing better detection of these metabolites in this study.

We also identified a new, previously undescribed potential ASD biomarker, homocitrulline. This metabolite was decreased in ASD patients and had the highest rank of all features in both SVM and PLS classification models. Homocitrulline is a poorly understood molecule which is known to be formed inside the mitochondria from lysine and carbamoyl phosphate. The decreased homocitrulline levels in the blood suggests that homocitrulline metabolism in the brain may also be disrupted, Homocitrulline levels are increased in urine and blood in patients with ornithine translocase (SLC25A15) deficiency which diverts carbamyl phosphate to react with lysine. These patients can exhibit behavioral abnormalities similar to ASD such as developmental delay, ataxia, spasticity, learning disabilities, cognitive deficits and/or unexplained seizures [57]. Rats treated with intracerebroventricular administration of homocitrulline are also observed to have disrupted brain redox status and energy metabolism [58], [59]. These observations suggest that elevated brain levels of homocitrulline are deleterious; however additional studies are needed to define the brain levels of homocitrulline and the potential role in the development of ASD.

Limitations

While the patient population was very well characterized, represented all severity levels of ASD and the blood samples were taken in a very systematic fashion, the total number of subjects (ASD = 52; TD = 30) was not large enough to definitively create a model for prediction of ASD or characterize metabolic differences that may be more highly associated with ASD subtypes.

Due to the small sample size, analysis of the data with respect to medication, sex, special diet, race, ethnicity or other potential confounding covariates was not conclusive. These are important considerations which require a larger sample size to address properly. A much larger study is currently in progress. Randomization methods were implemented to help prevent biasing the results based on any single covariate. The 4 misclassified patients in the 80 feature SVM model in the validation test set were Caucasian male and females whereas the Hispanic, Asian, and other ethnicity were predicted correctly.

Strong evidence of hemolysis, based on visual observation, was observed in 5 samples from ASD children. These samples were excluded from analysis Hemolysis can be a result of poor medical condition, but can also result from suboptimal blood collection and handling. However, Yin et. al. found that if EDTA tubes were placed on ice immediately, as was done in this protocol, the metabolome was very stable. Therefore, it remains possible, but unlikely, that minor hemolysis may confound the data by effecting some of the metabolite fold changes that were observed [60].

Future Endeavors

Further research is currently being carried out in much larger and younger patient populations to confirm these results, discover and confirm additional diagnostic metabolites and determine which are the most robust for evaluating ASD risk. We are also comparing metabolite profiles in clinically defined subtypes of ASD to determine whether predictive accuracy can be increased through better phenotyping of the ASD population. Analysis and consequential stratification will be performed based on covariables such as medication, sex, special diet, race, ethnicity, onset and co-morbid features such as gastrointestinal distress or seizure disorders may lead to a more accurate set of diagnostic metabolic profiles.