This retrospective study involved 1,253 cases who had GDD/ID and 100 typically developing children (TD) over the period from October 2015 to January 2021. Cases in the GDD/ID group were enrolled from the Xiangya Hospital, Central South University. Firstly, we screened positive results that could identify the cause of GDD/ID from 1,253 cases diagnosed with GDD/ID. Secondly, we screened the major differential metabolites and evaluated the diagnostic levels of differential metabolites in 132 urine organic acids from the remaining 1,230 negative results (863 GDD cases, 367 ID cases) that could not identify the cause of GDD/ID and 100 typically developing children (TD) (Supplementary Excel S1).

This study was approved by the Institutional Ethics Committee of Xiangya Hospital, Central South University. Both informed and written consents were obtained from the parents and/or legal guardians for study participation. Inclusion criteria for GDD/ID were as defined by the Diagnostic and Statistical Manual of Mental Disorders, Fifth Edition (DSM-V).

Exclusion criteria were: 1) presence of other diseases that may affect the development of the nervous system (e.g., intrauterine distress, history of hypoxia at birth, pathological jaundice, intracranial bleeding, perinatal infection, central nervous system infection, poisoning, and trauma history); 2) presence of chromosomal abnormalities, monogenic diseases, polygenic diseases which can cause GDD/ID. All participants were examined by senior pediatric neurologists. TD cases were enrolled from Children’s Health Clinic, Xiangya Hospital, Central South University.

Several precautions were strictly followed both before and after sampling to ensure specimen quality. The precautions and sampling steps were: 1) We took 10–20 ml of midstream urine in the morning, immersed 3 pieces of filter paper (the size was 5*5 cm) into the urine completely, took them out and dried them completely and put them in a specimen bag for inspection. Thereafter, we refrigerated it. Heating and drying were not allowed. 2) The qualified 3 specimens for each case were uniformly soaked and dried in urine filter paper sheets. 3) Samples were protected from the light and moisture and stored in refrigerator at 2–8°C for ≤ 14 days. 4) We rejected samples whose diaper infiltration area was too small, the specimen was contaminated or damaged, or moldy. The samples were collected from the inpatient or outpatient setting to ensure that external factors do not affect the samples. They were tested at the KingMed Diagnostics Laboratory (http://www.kingmed.com.cn/) (Guangzhou KingMed Diagnostics Group Co., Ltd, China). Gas chromatography/mass spectrometry (GC/MS) technology was performed for detection of organic acids. (See Supplementary Word S1 for specific test methods).

Information concerning study population characteristics was obtained from the Xiangya Hospital, Central South University electronic medical record system. Follow-up information was collected through regular clinic, WeChat and telephone communication. All assessments of children’s condition, daily life, and eating habits were provided by specialized clinicians and parents.

All urine organic acid test results were collected and entered into a spreadsheet database (Excel for Windows 10, Microsoft, China). Each item included identification number, name, sex, date of birth, test time, diagnosis, test report and results. R language script and IBM SPSS Statistics 22.0 software (IBM, Armonk, NY, United States) were used to analyse the data. Variables conforming to the normal distribution were expressed as mean ± standard deviation, and measurement data not conforming to the normal distribution was expressed as median or quaternary. Two-sample independent t-test was used for comparison. Categorical data were summarized in the form of frequencies and proportions, and analyzed with the Chi-square test or Fisher’s exact test where applicable. The level of significance was set at p ≤ 0.05.

Non-supervisory principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OLPS-DA) were used to compare the differences of major organic acids in urine between the affected children and TD group. We used a criteria of VIP>1.0 and p < 0.05 to filter out candidate metabolites. The differences of some organic acids in urine of each group were analyzed by heatmap. Metabo Analyst5.0, KEGG database, Metlin database and HMDB database were used to analyze relevant data and metabolic pathways. Receivers operating characteristic curve (ROC) and Youden index were used to determine the optimal diagnostic cut-off point of urine organic acids.

A total of 1330 participants were enrolled including 863 GDD cases (GDD group), 367 ID cases (ID group) and 100 TD children (TD group). A total of 1253 GDD/ID patients were enrolled including 23 patients with GDD/ID whose causes were determined by urine organic acid detection and 1230 patients with GDD/ID whose causes were unknown. The cause of 1.8% patients with GDD/ID could be identified by urine organic acid tests, including 7 cases of methylmalonate, 4 cases of glycerol kinase deficiency, 3 cases of glutaric aciduria, 2 cases of tyrosine metabolism disorder, 2 cases of Maple syrup urine disease. 2 cases of propionate aciduria, 1 case of ethyl malonate aciduria, 1 case of phenylketonuria, and 1 case of hyperlactic aciduria.

The urinary organic acids values of the cases in the GDD and ID groups were compared with those in the TD group, and the difference was statistically significant (p < 0.05). One hundred and thirty two kinds of organic acids were compared and the difference was statistically significant for 57 organic acids (Supplementary Table S1). Except for palmitic acid, the other 56 urinary organic acid values in GDD and ID groups were all higher than that of the TD group. The PCA results of the three groups with all acids are shown on Figures 1A,B. From the figure, we see that there are differences between the GDD group and the ID group and the TD group respectively. Most of the sample points fall within the 95% confidence interval of Hotelling’s T2 ellipse. (Supplementary Table S2). In order to eliminate the noise information that is not related to the classification, and also to obtain the related metabolite information that causes significant differences between the two groups, we used orthogonal partial least square discriminant analysis (OPLS-DA) to filter the signals that are not related to the model classification. The OPLS-DA results of the three groups are summarized in Supplementary Figures 1C,D and Supplementary Table S3. Both the GDD and the ID models had high explanatory and predictive rates, indicating that this model can distinguish between the two groups. OPLS-DA permutation plot results of GDD and ID groups showed no overfitting, and subsequent analysis of differential metabolite results was reliable. (Supplementary Figure S1)

According to the coefficient confidence interval of each variable and the VIP value of the variable weight, the meaningful difference variable (p < 0.05, VIP>1.0) was selected, that is, the specific biomarker. The identified potential marker metabolites are listed in Table 1. The comparisons of their levels between both ID and GDD groups versus TD group are also shown in Table 1. Compared with the TD group, the differential metabolites of the GDD group and the ID group had the same trend. Twenty four different metabolites were screened out in the GDD group, and 25 different metabolites were screened out in the ID group (VIP>1.0, p < 0.01). (Supplementary Excel S2) Except for palmitic acid, metabolites in GDD group and ID group were up-regulated, compared with TD group.

We calculated the Euclidean Distance Matrix for the quantitative values of the differential metabolites, clustered the differential metabolites by the complete linkage method, and displayed them in the form of heatmap (Figure 2). The outer circle of the circular heatmap represent different samples from GDD and TD groups, ID and TD groups, while the sides of the circular heatmap represent different metabolites. The color blocks at different positions represent the relative expression of metabolites at the corresponding positions. The redder the color, the higher the relative expression content, and the bluer the color, the lower the relative expression content. The top five identified possible biomarkers were palmitic acid, isocotonic acid, ethyl hydroxypropionic acid, ethylmalonic acid, and glycolic acid.

We calculated the correlation coefficient of the quantitative value of differential metabolites by Pearson method and presented it in the form of heatmap (Figure 3). We conducted correlation analysis on multiple differential metabolites to measure the degree of correlation between two differential metabolites, and the degree of correlation between two variables was represented by the correlation coefficient r. In the case of positive correlation, the value of r was between 0 and 1; in the case of negative correlation, the value of r was between −1 and 0. The closer the absolute value of r was to 1, the stronger the correlation between the two differential metabolites; the closer the absolute value of r was to 0, the weaker the correlation between the quantity variable.

The cut-off values of the different metabolites in the GDD and ID groups were calculated and sorted by using ROC curve and Youden index. Among the different metabolites in the GDD group, palmitic acid’s cut-off value was 7.190, area under the cut-off value curve (AUC) was 0.985, glycolic acid’s cut-off value was 2.210, AUC was 0.953, and 3-hydroxyisobutyric acid’s cut-off value was 2.315, AUC was 0.948. Among the different metabolites in the ID group, palmitic acid’s cut-off value was 5.690, AUC was 0.963, glycolic acid’s cut-off value was 2.810, AUC was 0.954, and 3-hydroxyisobutyric acid’s cut-off value was 2.330, AUC was 0.916. It can be seen from the Table 2 that the cut-off values of palmitic acid, glycolic acid, and 3-hydroxyisobutyric acid had high value for early identification of GDD/ID. The specificities of these three substances in the GDD group were 98.0, 90.0, and 91.0%, respectively, and for ID group were 99.0, 96.0, and 91.0%, respectively (Table 2). Supplementary Figure S2 shows the ROC curve of the cut-off value of palmitic acid, glycolic acid, 3-hydroxyisobutyric acid in the GDD group and the ID group. It can be seen from the figure that the AUC of all differential metabolites was higher than 0.9.

We tried to use bubble charts to show the possible related metabolic pathways (Figure 4). Each bubble in the bubble diagram represented a metabolic pathway. The abscissa of the bubble and the size of the bubble represented the size of the influence factor of the pathway in the topological analysis. The larger the size, the greater the influence factor. The vertical coordinate where the bubbles were located and the color of the bubbles represented the p value of enrichment analysis (take the negative natural logarithm, i.e. –In (P)). The darker the color, the smaller the p value and the more significant the enrichment degree. After bubble chart analysis, the pathway analysis of the differential metabolites in the GDD and ID groups were basically the same. According to the order of impact value and–In (P) value, the top 8 pathways were the synthesis and degradation of ketone bodies, citrate cycle (TCA cycle), alanine, aspartate and glutamate metabolism, pyrimidine metabolism, butanoate metabolism, pyruvate metabolism, fatty acid biosynthesis, valine, leucine and isoleucine degradation. This diverse distribution suggest that these organic acids may act on a variety of metabolic pathways and reflects the complexity of metabolic abnormalities in GDD/ID.

People constantly receive and respond to external or internal stimuli, and these experiences are learned and memorized in their brains. For this learning and memory process, nerve cells in the brain undergo enormous molecular and cellular changes, not only in the input -output -related local subcellular compartments but also in the central nucleus (Yu et al., 2011). The metabolic disorder of urine organic acid is caused by the lack of certain enzyme, which leads to the accumulation of related carboxylic acid and its metabolites, resulting in metabolic acidosis and damage of brain, liver, kidney, heart, bone marrow and other organs (Gabbi et al., 2017). Due to the accumulation of toxic metabolites, energy synthesis disorders, mitochondrial dysfunction, neuronal cell apoptosis, cytoskeletal phosphorylation changes and myelination disorders, patients' brain nerve structure damage, ganglioside and synaptic plasticity abnormalities, resulting in decreased learning and memory ability, and even intellectual disabilities (Gabbi et al., 2017).

The heatmap generated from our analysis of urinary organic acids showed the complex relationship among these compounds in GDD/ID. Several organic acids were in the same pathway, whereas others were involved in multiple pathways. To date, the metabolites that have been explored as possible GDD/ID biomarkers include: methylmalonic acid, propionic acid, methyl citrate, propionyl carnitine, propionyl glycine (Horster and Hoffmann, 2004; Wongkittichote et al., 2017), glutaric acid, 3-hydroxyglutaric acid (van Karnebeek and Stockler, 2012; Hoytema van Konijnenburg et al., 2021), 4-hydroxybutyric acid, γ-aminobutyric acid (Li et al., 2015), 3-hydroxyisobutyric acid (van Karnebeek and Stockler, 2012), 2-hydroxyglutaric acid (Kranendijk et al., 2012). The exact pathogenesis of the accumulation of these metabolites leading to GDD/ID is still unclear. The currently known pathogenesis mainly includes: energy metabolism disorders that affect structure and function of the brain (Flippo and Strack, 2017), energy substance supply disorders, neurotransmitter abnormalities that affect neurotransmission (Lim et al., 2020), the accumulation of harmful metabolic substrates or the inability of normal utilization of metabolites affecting the formation of myelin or lead to dysfunctional myelin (Marques and Saftig, 2019).

The OPLS-DA score plot showed a clear difference among the distribution characteristics of metabolites between GDD/ID cases and TD children. Our analyses showed that 5 urine organic acids had significant differences between GDD/ID cases and TD children and thus, could be valuable for future research on the etiology and pathogenesis of GDD/ID, and prognosticate diseases clinically.

The GDD group had higher levels of methyl citric acid, malic acid, 2-ketoglutaric acid, 2-ketone-isocaproic acid and 2-deoxy-4-hydroxyacetoacetic acid compared to TD children. Besides, the ID group had higher levels of methyl citric acid, malic acid, orotic acid, glycolic acid and 2-deoxy-4-hydroxyacetoacetic acid compared to TD children. On the contrary, both GDD group and ID group had decreased amounts of palmitic acid. These metabolites are associated with multiple biochemical processes (Koulman et al., 2009). The increase in orotic acid is more common in urea circulation disorders (UCDs), which leads to the dysfunction of the central nervous system (Savy et al., 2018). In our study, children with GDD/ID had a higher level of orotic acid in their urine, which was speculated to be related to UCDs, hyperammonemia, amino acid metabolism disorders, and ultimately cause damage to the nervous system. Palmitic acid is an intermediate product of fat acid synthesis and metabolism. Many key proteins related to neural development can be modified and processed by palmitoylation, which play important regulatory roles in the growth, polarization, and establishment of neural networks of neurons (Globa and Bamji, 2017). We speculate that reduction in palmitic acid content affects the palmitoylation processing level, thereby causing damage to the development of the nervous system.