The Children’s Autism Metabolome Project (CAMP, ClinicalTrials.gov Identifier: NCT02548442) study enrolled 1,102 children, ages 18–48 months, across 8 clinical sites from August, 2015 through January, 2018. We selected this age range because a consensus has emerged that a professional diagnosis of autism can be carried out accurately as early at 18 months of age. The centers included: The Children’s Hospital of Philadelphia; Cincinnati Children’s Hospital; The Lurie Center at Massachusetts General Hospital; The Melmed Center; The MIND Institute, University of California – Davis; Nationwide Children’s Hospital; The University of Arkansas for Medical Sciences; and Vanderbilt University Medical Center. Written informed consent from a parent or legal guardian was obtained and monetary compensation was provided to each participant. The study protocol was approved and monitored by Institutional Review Boards at each of the clinical centers.

Blood was collected from participants who had not eaten for at least 12 h by venipuncture into 6 mL sodium heparin tubes placed on wet ice (16, 26). Plasma was obtained after centrifugation (1,200 g for 10 min at room temperature) and stored at −80°C within 60 min of the blood draw. Hemolysis of samples was measured using spectrophotometry of the plasma (31). Plasma from hemolyzed samples with hemoglobin >600 mg/dL were excluded from metabolomics analyses and concentration values for xanthine, uric acid, or hypoxanthine were omitted when hemoglobin exceeded 300 mg/dL (26).

Three quantitative LC–MS/MS methods measuring 54 small molecule metabolites were performed in a CLIA-certified laboratory. The methods were analytically validated in compliance with FDA and CLSI guidance for bioanalytical method validation (32, 33). Quantification of analytes was performed using an Agilent Technologies G6490 triple quadrupole mass spectrometer. Detailed information about the sample preparation, detection, and quantification of metabolites can be found in the Supplemental Data Sheet. Analyte measurements below the lower limit of quantification (LLOQ) or above the upper limit of quantification (ULOQ) values were replaced with 90% of the LLOQ or 110% of the ULOQ value.

CAMP enrolled 1,102 participants and 916 met the inclusion and exclusion criteria described above (Supplementary Figure 1). Of these, 608 received a diagnosis of autism and 214 were considered TYP. The participant sample set was established after removing 32 autistic and 4 TYP samples that were hemolyzed, 77 autistic and 1 TYP participants’ samples that failed LC–MS/MS acquisitions, and 94 participants with developmental delay (DD) without autism. The final sample set contained 708 participant samples from 499 autistic and 209 TYP children.

We measured the concentrations of 54 metabolites and also evaluated the ratios of these metabolites (Supplementary Table 1). Metabolite ratio analysis can detect changes or reveal biological processes that may not be discerned by individual metabolites (34). For example, the concentrations of metabolites in a metabolic reaction sequence that has a minimal, but physiologically relevant, alteration of function of an enzyme or transporter may not show apparent alterations of the metabolites of that pathway. However, if the concentration of a metabolite that is distal to the metabolic reaction is expressed as a ratio to a metabolite that is proximal to the metabolic reaction, that ratio may reveal a change in that pathway that could otherwise go undetected. In addition, ratios of the concentration of metabolites can provide a normalization effect that reduces variation due to unrelated biological or analytical sources leading to higher specificity in diagnostic analyses (34). Clinical applications of metabolite ratio analysis include use in newborn screening for some inherited disorders of amino acid and organic acid metabolism and of mitochondrial fatty acid beta-oxidation (35, 36). Because of the benefits of ratio analysis, metabolite ratios are also utilized in association analyses with genetic variants and phenotypes to identify the metabolic underpinnings of more complex, multifactorial biological processes (37–39).

The metabolite ratios were generated from all unique combinations of the metabolites except 3-carboxy-4-methyl-5-propyl-2-furanpropanoic acid (CMPF), 4-ethylphenylsulfate (4-EPS) and dodecanedioic acid where 90% of the measurements were below the LLOQ. To create uniform visualization of metabolite ratios, the numerator and denominator were selected to yield a ratio with values that are greater than the biomarker threshold (see Biomarker Analysis). The concentrations of each metabolite or ratio of metabolites values were log base 2 transformed to reduce skewness and standardized (μ = 0, σ = 1) by subtracting the mean and dividing by the standard deviation prior to analyses. Participants’ metabolite measurements with missing data were omitted from analysis reducing the number of samples analyzed for a test statistic or imputed with the median value when missing values are not allowed by a test statistic. Analyses were conducted using R version 4.1.0 (40).

Receiver operator curve (ROC) analysis was used to select a biomarker value threshold (Figure 1) that maximized the percentage of autistic participants at a sensitivity above 4.5% when exceeded (26). The performance metrics were estimated using 4-fold cross-validation, repeated 50 times, stratified by participant sex, age, and diagnosis. This method of cross-validation trains and tests a model four times using independent sample sets, based on a training set of 75% and a test set of 25% of the samples, with model performance assessed as the average test set performance across repeats. Metabolites and ratios of metabolites were considered an autism-associated biomarker if the average performance had a sensitivity of at least 4.5% (indicating a subpopulation of at least 4.5% of the autism study participants that were biomarker-positive) and the proportion of the biomarker-positive (PMP) population of at least 90% autistic individuals (equivalent to the positive predictive value (PPV) of 90% within the CAMP study population prevalence). In addition to the sensitivity and PMP criteria, the permutation-based test statistic was significant at a false discovery rate adjusted p-value <0.1 (16, 26).

The final model thresholds were set using the entire study set of CAMP autism and TYP participants. These thresholds were used to generate participant biomarker outcome profiles by scoring a participant positive or biomarker negative for each biomarker. CAMP participant biomarker profiles were used to cluster the biomarkers, determine the prevalence of biomarker clusters in CAMP, and comparisons based on biomarker positive and negative populations. An overview of the methods used for biomarker selection and creating participant biomarker profiles is presented in Supplementary Figure 2. Optimization of biomarkers for identification of a likelihood of autism was performed using the process described previously (26).

Clustering was performed to reduce the complexity of the biomarkers by aggregating related biomarkers or biomarker clusters into groups. Hierarchical complete-linkage clustering of the participant biomarker outcome profiles was performed using the Jaccard distance based on scoring a biomarker outcome as 0 = negative and 1 = positive. The optimal number of biomarker clusters was estimated using the maximum value of the average silhouette width cluster validation index over a range of 5 to 50 clusters. The biomarker clusters were further evaluated by clustering the fold changes of metabolites by hierarchical complete-linkage clustering using a distance matrix based on the Pearson correlation coefficients (|1-r|) of the metabolite fold changes. The biomarker clusters dendrogram was cut based on similar patterns of fold changes for a subset of metabolites (see Results). Clustering was performed using the R packages ComplexHeatmap (41) and NbClust (42).

CAMP participant phenotypic and demographic information were based on physical and neurological examinations and behavioral testing performed by clinicians as well as a parental questionnaire. CAMP information from autistic and/or TYP children related to demographic information, diet, medications, behavioral assessments, or co-occurring conditions were selected for association analysis. The selected information was filtered to remove questions missing responses in more than 10% of autistic participants or that had an identical response in >98% of participants. The percent of ideal body weight (IBW) was based on the method of Traub and Johnson (43).

Association analysis of the biomarker values, biomarker defined subpopulations and the number of biomarker clusters to the CAMP demographic or phenotypic variables of autistic children was performed. To test for associations using biomarker values or the number of biomarker clusters, partial Spearman’s correlation coefficients rho (ρ) with age as a covariate were used when a demographic or phenotypic variable was continuous and a Kruskal-Wallis test was used when the metadata variable was categorical. Response wise Wilcoxon rank sum tests were used as post hoc tests for the Kruskal-Wallis tests. Association analysis between biomarker positive and negative populations and CAMP metadata was performed using Fisher exact tests for categorical variables and Welch t-tests for continuous metadata variables. Post hoc tests for Fisher exact tests of categorical metadata variables were performed by creating a dichotomous response variable for each response of the categorical variable followed by a response wise Fisher exact test. Effect sizes were reported using Cohen’s d for Welch t-tests, r statistic for Wilcoxon independent two sample tests, bias corrected Crammer’s V or odds ratios for Fisher’s Exact tests, and rank eta squared (η²) for Kruskal-Wallis tests. False discovery rate corrections (adj. p-value) were performed to control for multiple comparison testing (44). Partial correlations were calculated using the ppcor R package (45).

Plasma samples from 708 CAMP participants corresponding to 499 autistic and 209 typically developing children were utilized in this study. The autistic population was 2.5 months older than the TYP population (p-value <0.05, Table 1). The autistic population also had a higher proportion of males (79% versus 59%, Table 1). Asians and those who did not specify a race were overrepresented and Black or African Americans were underrepresented in the autistic population (Table 1). Other demographic factors found in Table 1 were balanced between the autistic and TYP populations. Parental interviews, medical records, and assessments by clinicians were used to evaluate co-occurring conditions associated with autism in the CAMP population as summarized in Supplementary Table 2.

We quantified 54 plasma metabolites comprised of organic and amino acids, acylcarnitines, purines and microbiome-associated metabolites that have been of interest in autism research (15–17, 20, 26, 46, 47) (Supplementary Table 1). We also evaluated 1,275 unique ratios of these metabolites. Our aim was to test if a quantitative threshold of the concentration of a metabolite or the value of a ratio of metabolites could serve as a biomarker of an autism-associated subpopulation (Figure 1). Of the metabolites and ratios of metabolites that were examined, 153 autism-associated biomarkers were identified. The biomarkers included 4 individual metabolites and 149 ratios of metabolites (Supplementary Table 3). Each biomarker identified a subpopulation that contained between 4.5% and 11% of the CAMP study autistic population. The subpopulations were comprised of 90–100% autistic participants. One or more biomarkers identified 414 (83%) of the CAMP study autistic participants. The most frequent metabolites to appear as biomarker numerators or denominators were lactate, carnitine, pyruvate, leucine, glycine, octanoylcarnitine, citrate, 4-hydroxyproline, phenylalanine, and 2-ketoglutarate. The metabolites 3-aminoisobutyric acid, serotonin, inosine, 3-indoxyl sulfate, 4-ethylphenyl sulfate, CMPF, dodecanedioic acid, hydroxybutyrylcarnitine, indoleacetic acid, indolelactic acid and p-cresol sulfate were not identified as biomarkers.

An autistic participant is often identified by more than one of the 153 biomarkers (median 6 biomarkers, range 1–62, Figure 2, columns). The identification of an individual by multiple biomarkers suggested that the biomarkers may be identifying related metabolic processes (Supplementary Table 4). To explore this hypothesis and reduce complexity of the biomarkers, the biomarkers that often co-identify participants were grouped together using hierarchical clustering (Figure 2 rows). The 153 biomarkers formed 30 clusters; 23 contained 2 to 40 biomarkers and 7 clusters each consisted of 1 biomarker (Figure 2 rows, Supplementary Table 5). The biomarker clusters ↑Glycine, ↑Ornithine, ↑Gly↓BCAA, ↓Leucine, ↑LacPyr, ↑LacPyr↓GlyAsp, ↓aKG1, ↓Urate, and Arg/4hyp contain additional metabolites and ratios of metabolites that were not evaluated in our earlier work (16, 26). New biomarker clusters include those that contain the metabolites malate, citrate, xanthine, hypoxanthine, taurine, carnitine and the indicated acylcarnitines.

The autism-associated biomarkers can either be used by themselves or clustered with other autism-associated biomarkers to identify autism-associated metabolic phenotypes. The 30 biomarker clusters each identified subpopulations with a range of 5% to 28% of the CAMP autistic participants and were each comprised of 89% to 97% of autistic individuals (Table 2). As with individual biomarkers, participants are often identified by more than one biomarker cluster (median 3 biomarker clusters, range 1–19) with 52% of study participants identified by more than one biomarker cluster. The observation that different biomarker clusters often co-identify an individual (Supplementary Figure 3; Supplementary Table 6) suggests that the metabolic clusters are not independent of each other, i.e., there are alterations of metabolism that commonly occur together in the subpopulations.

To investigate the relationship of biomarker clusters and metabolism, we compared the fold changes of metabolites between biomarker cluster-positive and negative participants. We evaluated the changes in metabolites to obtain a more precise understanding of the metabolic processes and pathways associated with biomarker clusters that could not be derived from ratios of metabolites. Metabolites with the greatest differences in concentrations between cluster positive and negative populations provide information about the influence of a metabolite in a biomarker cluster’s metabolic phenotype. Patterns of metabolite changes among the biomarker clusters can then be used to group biomarker clusters based on the similarity of fold changes across metabolites. We grouped the biomarker clusters fold changes and selected 7 groups of biomarker clusters based on acylcarnitines, branched chain amino acids, lactate, and pyruvate (Figure 3; Supplementary Table 7). By focusing on metabolites with the greatest changes and grouping the biomarker clusters, we were able to reduce the complexity of our analysis and obtain more meaningful biochemical interpretation of the biomarker clusters which are described below.

Four groups of biomarker clusters had different patterns of fold changes where increased acylcarnitines were evident among them (Figure 3). Group 1 consists of the ↑LacPyr, aKG/Val, ↑LacPyr↓GlyAsp, and ↑Malate clusters and is present in 35% of autistic participants. The group is characterized by increases in lactate, pyruvate, and malate, with varied increases in saturated medium and long chain acylcarnitines. A second pattern of fold changes is evident in the ↑C10C8C6 cluster of group 2 and is observed in 19% of the autistic participants. This cluster’s fold changes are characterized by strongly increased medium chain acylcarnitines, elevated palmitoylcarnitine and decreased glutamate and aspartate. The third pattern of fold changes with increased acylcarnitines consists of two clusters, Leu/Car and Leu/Val, is present in 11% of participants with autism and is characterized by increased levels of hexanoylcarnitine, palmitoylcarnitine and the branched chain amino acids (BCAAs) and decreased carnitine.

Two groups of biomarker clusters had a pattern of fold changes with decreased levels of acylcarnitines (Figure 3, groups 5 and 7). The clusters ↓Carnitine1, ↑Glycine, ↓Glutamine, ↑Arginine, ↑Gly↓BCAA, ↓Leucine, that together comprise 37% of the autistic participants, are characterized by decreased levels of the BCAAs, medium and long chain acylcarnitines and glutarylcarnitine, in combination with increased alanine, glycine, and ornithine. Another pattern of fold changes is present in 50% of autistic participants, shared by the ↓Urate, Lys/Gln, ↓aKG2, ↑Alanine, ↑LacPyr↓(C6-C16), ↓C5DC, ↑Ornithine, Lys/Phe, ↓aKG1, and ↓Citrate clusters, had decreased medium and long chain acylcarnitines, decreased glutarylcarnitine, variably decreased TCA cycle metabolites and variably increased alanine, glycine, threonine, methionine, proline and ornithine.

The biomarker clusters ↓Hydroxyproline, ↓Carnitine2, ↑Taurine, Eta/Kyn, ↑TauXan↓C16, and ↑Ethanolamine in group 4 (Figure 3, group 4) exhibit a pattern of fold changes occurring in 40% of the autistic participants. Increases of aspartate, glutamate, hypoxanthine, xanthine, taurine, ethanolamine, lactate and pyruvate are associated with this cluster. A pattern of fold changes noted in a single cluster, Arg/4Hyp (Figure 3, group 6), shows decreased 4-hydroxyproline and tricarboxylic acid cycle (TCA) intermediates and increased arginine. These clusters do not exhibit consistent changes in the acylcarnitines.

In terms of specific metabolites, decreased free carnitine is noted in three clusters that appear in different groups, ↓Carnitine1, ↓Carnitine2 and Leu/Car (Figure 3, groups 3–5). These biomarker clusters are observed in 26% of autistic participants. Each of these is associated with a distinctive metabolic profile. The first has decreased medium and long chain acylcarnitines and increased alanine, glycine, ornithine and proline. The second has increased lactate, pyruvate, malate and hypoxanthine and minimal changes of acylcarnitines. The third is characterized by increased BCAAs, hexanoylcarnitine and palmitoylcarnitine.

The biomarkers use quantitative thresholds to identify potentially diagnostic subpopulations of the autistic participants that are comprised of at least 90% autistic individuals. The thresholds provide an opportunity for the biomarkers to be used as “tests” for biomarker levels that are associated with a likelihood of autism. The tests individually have high specificities (>98%) but low sensitivities (5–11%). The high specificity of these tests allows for a ‘stacking approach’ whereby a subset of the 153 biomarker-based tests can be combined into a battery of tests that increases overall test sensitivity for a likelihood of autism. Using a test optimization process to maintain a specificity of at least 90% (26), a subset of 42 biomarkers was selected that identified CAMP autistic participants with a sensitivity of 72% and specificity of 90% (Supplementary Table 8).

The levels of biomarkers corresponding to metabolite concentrations and ratios of metabolite values from autistic participants were evaluated for associations with behavioral test scores, demographic factors, co-occurring conditions, diets, and medications.

The biomarker-positive autistic individuals identified by the individual biomarkers or biomarker clusters were tested for associations with CAMP autistic participants’ metadata. Participants identified by the biomarkers Methionine/Carnitine, Ornithine/Carnitine, Alanine/Tyrosine, Alanine/Phenylalanine or Ethanolamine/4-hydroxyproline were associated with lower DQ and MSEL subdomain scores (Supplementary Table 11). Individuals identified by Leucine and the biomarker clusters ↑LacPyr or ↑Malate were associated with a decreased MSEL Receptive Language Score, but not a decrease in DQ (Supplementary Table 11). Participants identified by Taurine/Malate were associated with lower BMIs (Supplementary Table 11). Other metadata were not significantly associated with participants identified by an individual biomarker or a cluster of biomarkers at a significance threshold of an adjusted p-value of 0.1 (Supplementary Table 12).

The autistic participants were evaluated for associations with autism severity based on the ADOS-2 test scores (27) and the number of biomarker clusters that identify an autistic participant (Figures 4A–C; Supplementary Table 13). Increased ADOS-2 CSS, indicating more severe autism characteristics, was positively associated with the number of biomarker clusters (ρ = 0.14, adj. p-value = 0.016). The number of biomarker clusters was also correlated (ρ = 0.12, adj. p-value = 0.05) with more severe SA CSS. The number of clusters was not associated with RRB CSS (ρ = 0.05, adjusted p-value = 0.32). That the number of clusters have an association with SA CSS and lack of association with RBB CSS suggests that the increased autism severity associated with increased number of clusters is related to deficits in social communication.

The autistic participants with increasing numbers of biomarker clusters exhibited more severe developmental delays as measured by lower MSEL test scores (Figures 4D–H; Supplementary Table 13). The DQ decreased with increasing numbers of clusters (ρ = − 0.13, adj. p-value = 0.008). Decreased DQ was associated with a negative correlation of expressive (ρ = −0.14, adj. p-value = 0.005) and receptive language (ρ = −0.14, adj. p-value = 0.005), fine motor (ρ = −0.13, adj. p-value = 0.013) and visual reception (ρ = − 0.12, adj. p-value = 0.016) with increasing numbers of clusters.

The number of biomarker clusters was not associated (Kruskal-Wallis test adjusted p-value > 0.1 or Spearman partial correlation test adjusted p-value>0.1) with age, race, BMI, diet, medications, sleep problems, gastrointestinal issues, food selectivity, or delayed attainment of early language or motor milestones (Supplementary Table 14). The lack of association with these factors raises the interesting possibility that an increased number of biomarker clusters is primarily associated with the severity of autism and reduced developmental quotients in the CAMP autistic population.