Changes in the Non-targeted Metabolomic Profile of Three-year-old Toddlers with Elevated Exposure to Polycyclic Aromatic Hydrocarbons

Polycyclic aromatic hydrocarbons (PAHs), which are widely distributed in the environment owing to the incomplete combustion of fossil fuels, organic matter pyrolysis, and geochemical processes, have recently gained considerable attention. In addition to their well-established carcinogenic and cardiovascular toxicities, there has been increasing focus on the neurodevelopmental toxicity associated with PAH exposure. Evidence suggests a positive association between high levels of PAH exposure during early childhood (5–7 years) and symptoms related to anxiety/depression as well as attention problems^[1,2]. According to the World Health Organization, neurobehavioral disorders, such as depression and autism, have emerged as the third most prevalent diseases after cancer and cardiovascular diseases, and depression-related suicide is the second leading cause of death among adolescents^[3]. Therefore, addressing the neurodevelopmental damage caused by PAHs has become an imperative public health concern that cannot be overlooked.

Animal studies suggest that exposure to PAHs during pregnancy may result in delayed neurodevelopment, lower IQ, abnormal behavior, and other adverse health effects. This association is also supported by epidemiological studies^[1,4]. Importantly, the human brain undergoes rapid development not only during the third trimester of pregnancy but also until at least two years after birth^[5]. Our previous cross-sectional investigation^[6] and cohort study further indicated that PAH exposure within the first two years after birth could contribute to autism-like neurodevelopmental and behavioral changes^[7].

Along with rapid brain development, the first two or three years of early life play a crucial role in establishing and maturing the gut microbiota in infants’ intestines, resembling that of adult flora^[8]. Humans are exposed to PAHs via inhalation, gastrointestinal ingestion, and dermal absorption. Young children, particularly infants and toddlers, are primarily exposed to PAHs through diet^[9] and unintentional dust ingestion. Notably, unintentional dust ingestion resulting from hand-to-mouth behavior can contribute significantly (up to 42%) to non-dietary PAH exposure^[10]. Owing to their high lipid solubility, PAHs are easily absorbed into the gastrointestinal tract of mammals^[11]. Therefore, PAHs can easily reach and affect the gastrointestinal tract. Persistent organic pollutants (POPs) affect the diversity of the intestinal microbiota by activating aromatic hydrocarbon receptors (AhRs), thereby affecting host homeostasis^[12] and contributing to the regulation of behavioral and physical abnormalities associated with neurodevelopmental disorders. As a significant category of POPs, PAHs serve as important ligands for AhRs. Our pilot study demonstrated that PAH exposure was negatively correlated with most Bacteroidetes species but positively correlated with Firmicutes abundance, which is consistent with other animal studies^[13]. Additionally, in vitro observations of human fecal microbiota exposed to benzo[a]pyrene (BaP) revealed dose-dependent alterations in gut bacterial activity related to the volatile metabolome and transcriptome. The activated pathways were associated with vitamin and aromatic compound metabolism, cell wall chemical synthesis, and DNA replication and repair. The inhibited pathways included microbial chemotaxis and hydrocarbon and glucose metabolism^[14]. Another type of PAH, phenathrene, was found to increase the Firmicutes/Bacterioidetes ratio, leading to hepatic lipid accumulation. However, hepatic lipid metabolism was restored after fecal microbial transplantation. These findings emphasize the importance of establishing a correlation between the composition or functional activity of the gut microbiota and changes in the host metabolic profiles. Nevertheless, to the best of our knowledge, there is currently no study on the alterations in urinary metabolic profiles caused by PAH exposure in toddlers.

Owing to the variety of sources and types of PAHs, relying solely on a single biomarker, such as 1-hydroxypyrene (1-OHPy), is insufficient for evaluating the total PAH exposure. To comprehensively assess individual PAH exposure, researchers have attempted to detect multiple hydroxyl metabolites, including 1-hydroxynaphthalene (1-OHNap) and 2-hydroxynaphthalene (2-OHNap) from naphthalene, 2-hydroxyfluorene (2-OHFlu) from fluorine, 1-, 2-, 3-, 4-, 9-hydroxyphenanthrene (1-,2-,3-,4-,9-OHPh) from phenanthrene, and 1-OHPy from pyrene^[15]. Metabolomics has the potential to elucidate lifestyle, environmental, and genetic factors^[16], facilitating the identification of new intermediate biomarkers and their application in molecular epidemiology, owing to its comprehensive information content^[17]. The combination of ultrahigh-performance liquid chromatography and quadrupole time-of-flight mass spectrometry (QTOF MS) is widely used in metabolomic research. Hydrophilic interaction chromatography (HILIC) is recognized for its complementary selectivity to reversed-phase chromatography. Consequently, we investigated whether PAH exposure could alter the metabolite profile of noninvasive urine samples from 3-year-old toddlers using HILIC-ESI (±)-QTOF MS. Additionally, we explored the potential of these identified metabolites as biomarkers for PAH exposure and their correlation with gut microbiota implicated in neurodevelopment.

OxiSelect™ BPDE DNA Adduct ELISA Kit was purchased from Cell Biolabs (San Diego, Canada). Ten hydroxyl metabolites of PAHs (1-OHNap, 2-OHNap, 2-OHFlu, 1-, 2-, 3-, 4-, 9-OHPh, 1-OHPy, and 6-hydroxychrysene), mass spectrometry grade, were purchased from Sigma-Aldrich (St. Louis, Missouri). QIAamp DNA Stool Mini Kit was purchased from Qiagen (Germantown, MD). Acetonitrile and methanol (HPLC grade) were purchased from Merck (Darmstadt, Germany).

Participants were randomly selected from a birth cohort established in 2014 that included children born between January and December of that year at various hospitals in Qingdao^[18]. A total of 1,973 children were included in this study. Subsequently, children up to the age of three within this cohort were screened for eligibility, excluding those with certain medical conditions or recent use of antibiotic/microecological regulators. The exclusion criteria included cardiovascular diseases, metabolic diseases, obesity, inflammatory bowel diseases, irritable bowel syndrome, asthma, previous organic brain diseases, and special diets. Additionally, those who had taken antibiotics or microecological regulators within two months prior to sampling were excluded. Among the initial cohort population with cord blood BaP-DNA adduct and urine OH-PAH determination at three years old (n = 400), 20 children with the highest total concentration of the targeted OH-PAHs and 20 children with the lowest concentration were selected as representatives of the high- and low-exposure groups, respectively. After excluding three individuals with abnormally high concentrations of BaP-DNA conjugate, 18 children remained in the high-exposure group and 19 remained in the low-exposure group for the final analysis. The OH-PAHs levels in the two groups are presented in Supplementary Table S1 (available in www.besjournal.com). This study was approved by the Ethics Committee and the Institutional Review Board of Qingdao Center for Disease Control and Prevention (No. 201802). Written informed consent was obtained from legal guardians of all participants involved in this research project. The authors declare that all procedures complied with ethical standards and adhered to the guidelines of Helsinki Declaration of 1975, as revised in 2008.

We conducted a follow-up questionnaire survey on the established birth cohort. The questionnaire consisted of two sections: maternal and neonatal situation questionnaires. The maternal questionnaire covered age, height, weight, educational level, health history, occupation, living environment, dietary habits, smoking habits, and alcohol consumption patterns during pregnancy. The neonatal questionnaire included gestational week, birth weight, mode of delivery, incidence of neonatal asphyxia, and occurrence of birth injury. Two months before sampling, a supplementary questionnaire survey was conducted to eliminate potential confounding factors that could affect the intestinal flora. The survey focused on cardiovascular diseases, metabolic diseases, obesity status, inflammatory bowel diseases, irritable bowel syndrome, asthma, previous organic brain diseases, special diet, and recent use of antibiotics or microecological regulators.

The neurobehavior of children aged 2–3 years was assessed by professional child healthcare doctors using the Achenbach Child Behavior Checklist (CBCL)^[19]. The CBCL consists of 100 items that evaluate neurobehavioral abnormalities in six behavioral symptom factors: anxiety/depression, withdrawal, sleep problems, physical problems, aggressive behavior, and destructive behavior. Internalizing behavior was defined as the sum of anxiety/depression and withdrawal factors, whereas externalizing behavior was defined as the sum of aggressive and destructive factors.

The intellectual development assessment for children was conducted by a single trained medical professional using the Gesell Development Diagnosis Scale to ensure consistency across scales. The assessment covered five behavioral domains: adaptive, gross motor, fine motor, language, and personal-social behavior. The test results are expressed as the development quotient (DQ)^[20].

The collection of umbilical cord blood and the analysis of PAH-DNA adducts have been described in detail in our previous studies^[7].

Morning first-void urine samples (15 mL) were collected consecutively for 5 days from 3-year-old children, refrigerated, and sent to the designated hospital by their parents. After thorough mixing, the samples were frozen and stored at −70 °C in the Qingdao Center for Disease Control and Prevention for subsequent analysis. Urinary hydroxyl metabolites of PAH (OH-PAHs) were tested following the method described by Guo et al.^[15] with minor modifications^[18]. The OH-PAHs selected in this study included ten exposure biomarkers: 1-, 2-OHNap; 2-OHFlu; 1-, 2-, 3-, 4-, and 9-OHPh; 1-OHPy; and 6-OHChr (Supplementary Materials, available in www.besjournal.com).

After slowly thawing the morning first-void urine sample at 4 °C, it was mixed with 100 μL of pre-cooled methanol/acetonitrile solution (1:1, v/v) and vortexed. The mixture was then left to stand at −20 °C for 60 min. Subsequently, centrifugation was performed at 14,000 ×g and 4 °C for 20 min to obtain the supernatant, which was collected and subjected to vacuum drying. For mass spectrometry analysis, the dried sample was reconstituted by adding 100 μL of acetonitrile aqueous solution (acetonitrile: water = 1:1, v/v), followed by vortexing and centrifugation (14,000 ×g, 4 °C, 15 min) to collect the supernatant. To ensure stability and repeatability in instrument analysis monitoring, quality control (QC) samples were prepared by pooling together a volume of each individual sample (10 μL); these QC samples were regularly inserted into every eighth set of analyzed samples.

Urine samples were analyzed using an Agilent 1,290 Infinity LC system (Agilent Technologies, Santa Clara, California) coupled with an AB SCIEX Triple TOF 5600 System (AB SCIEX, Framingham, MA) at Shanghai Applied Protein Technology Co., Ltd. (http://www.aptbiotech.com). Chromatographic separation was performed on ACQUITY UPLC BEH Amide 1.7 µm (2.1 mm × 100.0 mm) columns (Waters, Co., Ltd.) in both positive and negative modes. Detailed chromatographic and mass spectrum conditions can be found in the Supplementary Materials section (Supplementary Materials). Spectrogram comparison and Principal component analysis (PCA) were used for QC and assurance of the QC samples (Supplementary Materials).

Fresh fecal samples were collected at three years of age and transported to the Qingdao CDC laboratory with ice packs within 2 h of collection. All samples were immediately frozen and stored at −80 °C before analysis. DNA from the gut microbiota was extracted from the stool samples using the QIAamp DNA Stool Mini Kit. Subsequently, the DNA samples were sent to the TinyGene Company (Shanghai, China) for 16S rRNA gene sequencing analysis using Illumina MiSeq (V4-V5 region).

IBM SPSS 21 (IBM Inc., Chicago, IL) was used for data analysis. The CBCL scale utilizes the actual score as a measure of children’s behavior. Metabolites with over 50% missing values were excluded from the dataset and not included in subsequent analyses. The means of the total hydroxyl PAH metabolites between the two groups were compared using Student’s t-test, while the Mann–Whitney test was employed to compare the CBCL and DQ scores of the two groups. The total peak area of the metabolic data was normalized and Pareto-scaled using SIMCA-P software. PCA was then performed on the peaks extracted from all samples and QC samples after Pareto scaling (Supplementary Figures S1-S2, available in www.besjournal.com). The closely clustered QC samples in both the positive and negative ion modes indicated excellent repeatability of the experimental procedures conducted in this project. The OPLS-DA method, a supervised discriminant statistical approach that is well-suited for distinguishing differences between two groups and has become the standard procedure in metabolomics^[21] was used to evaluate the differences between the two groups. Based on the variable importance for projection (VIP) obtained using the OPLS-DA model, we measured how the expression pattern of each metabolite influenced the classification and discrimination of the group samples, allowing us to determine differential metabolites with biological significance. In this study, VIP > 1 and P < 0.1 were used as criteria for screening out differential metabolites, where P < 0.05 represented a significant difference. Furthermore, a permutation test was conducted to assess the robustness of the model. Metabolites with significant differences were subjected to hierarchical clustering and KEGG analysis (https://www.genome.jp/kegg/). The correlation analysis results were expressed as r-values. A two-tailed test with P < 0.05 indicated statistical significance. To gain insights into the data structure, clustering heat maps and networks were generated using the R software.

Figure S1.  TIC spectrogram overlap comparison. The above was for positive ion mode and below was for negative ion mode.

Figure S2.  PCA score plot for metabolomics analysis of urine samples from the high and low PAHs exposure groups and QC samples. (A) Positive ion mode; (B) Negative ion mode.

The high- and low-PAH exposure groups did not differ significantly in terms of sex, maternal age, gestational age, or birth weight. However, a significant discrepancy was observed in the urinary OH-PAH concentrations (P < 0.001) between the two groups (Table 1, Supplementary Table S1). No variation was found in the BaP-DNA adduct levels in the umbilical cord blood samples. Additionally, no notable differences (P > 0.05) were identified through neurodevelopmental tests for either DQ or CBCL syndrome scores between the two groups (Table 1).

In addition to the significant difference in total OH-PAHs between the high- and low-exposure groups (30.63 vs. 3.83; P < 0.001) (Table 1), the hydroxyl metabolites of naphthalene, fluorene, and phenanthrene were also significantly higher in the high-exposure group than in the low-exposure group (P < 0.01) (Supplementary Table S1).

The OPLS-DA model parameters obtained through 7-fold cross-validation are shown in Figure 1. A permutation test was conducted to assess the robustness of the model, and no evidence of overfitting was observed (R² = 0.759–0.779, Q² < 0.5) (Supplementary Figure S3, available in www.besjournal.com). The OPLS-DA model effectively distinguished between the samples with high and low PAH exposure (Figure 1). Although the interpretability of the model (R²Y > 0.5) was satisfactory, its predictability (Q² < 0.5) was poor.

Figure 1.  OPLS-DA score plot for metabolomic analysis of urine samples from the high and low PAH exposure groups. (A) Positive ion mode. (B) Negative ion mode.

Table 2 presents 15 metabolites with significant differences between the positive and negative ion modes. Figure 2 displays the clustering analysis results of differential metabolites in the high- and low-exposure groups under positive and negative ion modes. Table 3 shows that the hydroxyl metabolites of naphthalene, fluorine, and phenanthrene were significantly associated with some differential metabolites (r = 0.336–0.531). Although no significant changes in the antioxidant uric acid level were observed, there were increased concentrations of upstream metabolites, including adenosine (VIP = 1.38, P = 0.063), deoxyinosine (VIP = 1.56, P = 0.089), xanthine (VIP = 1.46, P = 0.043), and hypoxanthine (VIP = 4.05, P = 0.077). Conversely, the concentration of downstream urea decreased (VIP = 1.00, P = 0.096). Three other modified nucleotides, N4-acetylcytidine, 2-methylguanosine, and pseudouridine, were also elevated (Figure 2). Five amino acid metabolism-associated metabolites (D-glucuronate, L-gulonic gamma-lactone, L-threonate, hippuric acid, and gentisic acid) exhibited significant differences between the groups (P < 0.05) (Table 2, Figure 2). The concentration of glycerophosphocholine was positively associated with PAH exposure (P < 0.01) (Table 2, Figure 2), particularly with the metabolites of naphthalene (r = 0.377) and phenanthrene (r = 0.485) (Table 3). We also observed that the pantothenate level was significantly higher in the high PAH exposure group than in the low PAH exposure group (P < 0.01) (Table 2, Figure 2), and it exhibited significant associations with the metabolites of naphthalene (r = 0.386), fluorine (r = 0.474) and phenanthrene (r = 0.380) (Table 3).

Figure 2.  Clustering analysis of differentially expressed metabolites in the high- and low-exposure groups. (A) Positive ion mode; (B) Negative ion mode.

Differentially expressed metabolites were analyzed using KEGG (https://www.genome.jp/kegg/) to determine their relevant pathways and biological functions. The majority of differential metabolites were concentrated in amino acid (such as aromatic amino acids [AAAs], arginine, proline, glycine, lysine, alanine, aspartic acid, serine, glutamic acid, taurine, and hypotaurine), carbohydrate (such as amino sugar, nucleoside sugar, inositol phosphate, and pentose phosphate), lipid (such as glycerophospholipid, sphingolipid), and nucleotide metabolic pathways (such as purine and pyrimidine) and vitamin-related processes (such as biotin, folic acid, pantothenate, and vitamin B6). Additionally, differential metabolites also participated in certain signal transduction pathways (such as cAMP and cGMP-PKG signaling pathways) involved in environmental and genetic information processing, and played a role in disease-specific metabolism, such as neurodegenerative diseases (Parkinson’s disease) (Supplementary Table S2, Supplementary Figure S4, available in www.besjournal.com).

The correlation network, based on Spearman’s correlation coefficients, revealed a positive association between increased levels of deoxyinosine, hypoxanthine, and N4-acetylcytidine and the abundance of the Firmicutes genera Subdoligranulum and Coprococcus (r = 0.36–0.45) (Figure 3). However, no significant correlation was observed at the phylum level (Supplementary Figure S5, available in www.besjournal.com). Similarly, the concentration of 2-methylguanosine was significantly positively correlated with the abundance of Subdoligranulum (r = 0.50) and Coprococcus (r = 0.44) but was negatively correlated with the genus Bacteroidetes from the Bacteroides phylum (r = –0.45). The levels of both pantothenate and glycerophosphocholine were positively correlated with the anti-inflammatory bacterial taxon Faecalibacterium, a genus belonging to the Firmicutes phylum, with r values of 0.46 and 0.38, respectively. Pantothenate was also positively correlated with certain bacterial taxa from the Bacteroidetes and Proteobacteria phyla, with an r value of 0.35. L-gulonic gamma-lactone was positively correlated with the Fusobacterium genus of Fusobacteria (r = 0.45) but was negatively correlated with the Flavonifractor genus of Firmicutes (r = –0.45). A positive association was observed between Coriobacteriales of the Actinobacteria phylum and acetyl phosphate (r = 0.40–0.45). Most of the correlations between metabolites and genera remained consistent across taxonomic levels (Supplementary Figures S5–S9, Supplementary Tables S3–S4, available in www.besjournal.com).

Figure 3.  Microbiota (genus)-metabolites correlation network based on Spearman’s correlation coefficients. Each node represents one genus (circle) or metabolite (triangle), and two nodes are linked if the correlation is significant (two-sided pseudo P < 0.05) and > 0.3. Lines between nodes indicate positive correlations (red lines) or negative correlations (blue lines). The node size is proportional to the number of correlated metabolites. (A) positive ion mode; (B) negative ion mode.

This study suggests that high PAH exposure in young children may alter their urinary metabolomic profiles, affecting various metabolic pathways such as those of amino acids, lipids, carbohydrates, nucleotides, and vitamins. However, these differential metabolites were only partially effective in predicting individual exposure. Additionally, a significant number of these altered metabolites were associated with inflammatory responses and were potentially influenced by the intestinal flora, which could impact neurodevelopment.

Our study suggests that early exposure to PAHs in young children can lead to changes in systemic metabolic homeostasis, although PAH exposure could not be predicted using non-targeted metabolites in children’s urine. Similarly, a systematic review found strong correlations at the group level between PAH metabolites and PAH-DNA adducts in relation to environmental exposure to PAHs, even at low levels of air pollution caused by BaP; however, these biomarkers have limited utility when assessing individual-level environmental exposure^[22]. Therefore, although our study successfully distinguished between high and low PAH exposure using non-targeted metabolites, these metabolites were not predictive of individual PAH levels.

The involvement of systemic inflammation and oxidative stress in cardiovascular diseases^[23], respiratory diseases^[24], and metabolic syndromes^[25] caused by PAHs has been reported even in healthy individuals exposed to PAHs^[26]. Urinary hydroxyl metabolites of phenanthrene and fluorene were positively associated with serum C-reactive protein levels in the National Health and Nutrition Examination Survey^[27], linking PAH exposure to systemic inflammation and oxidative stress. Additionally, evidence suggests that PAH exposure is involved in inflammation or oxidative stress in specific organs such as the gastrointestinal tract^[28]. The destruction of cellular signaling pathways associated with oxidative stress or inflammation serves as the primary driving force behind alterations in metabolic pathways involving various compounds such as nucleotides, amino acids, lipids, and bile acids^[29].

The results showed that PAH exposure induced abnormal purine metabolism. A previous study showed that low levels of PAH exposure in the general population can cause oxidative stress, resulting in abnormalities in the TCA cycle and amino acid metabolic pathways, leading to mitochondrial dysfunction and inadequate energy supply; simultaneously, oxidative stress or inflammatory responses may also trigger abnormal lipid metabolism^[30]. Therefore, it appears that PAH exposure induces systemic inflammation or oxidative stress, which leads to a decrease in antioxidant uric acid concentrations and its downstream product, urea, while concurrently increasing the concentrations of upstream metabolites such as adenosine, deoxyinosine, oxanthine, guanosine, and hypoxanthine. Our findings are consistent with those of a previous study on BaP exposure in lung cancer H460 cells, where BaP levels were similar to those observed in human blood at an equivalent concentration range (5 µg/L). BaP-exposed lung cancer cells showed a significant accumulation of purine nucleotides, including adenosine monophosphate, adenosine 5′-diphosphate, inosine, hypoxanthine, guanine, and guanosine. This suggests an altered metabolism and/or transport of purine nucleotides following BaP exposure^[31].

Elevated levels of other modified nucleotides, such as N4-acetylcytidine, 2-methylguanosine, and pseudouridine, may also be associated with inflammatory responses. The urine concentration of N4-acetylcytidine was reported to increase by 24% in patients with interstitial cystitis and was inversely correlated with uromodulin levels^[32]. Uromodulin has been implicated as a potential trigger for the activation of inflammatory cytokine interleukin-1β (IL1β)-dependent innate immunity^[33]. Additionally, the plasma levels of N4-acetylcytidine were significantly lower than the standard levels in 25 patients with pulmonary fibrosis^[34]. Furthermore, N4-acetylcytidine can stimulate neuroglia and induce high-mobility group box 1 signaling to sustain an NLRP3 neurogenic inflammatory response^[35]. Increased levels of 2-dimethylguanosine and pseudouridine have been observed during oxidative stress in patients with pneumonia, sepsis, cancer, and cardiovascular diseases^[36]. Furthermore, under conditions of inflammatory stress, enzymatic post-transcriptional modification of RNA can occur, leading to increased release of modified nucleosides^[37]. Pseudouridine is abundantly present in both rRNA and tRNA, exerting its role as a damage-associated molecule by stimulating the expression of inflammatory cytokines such as TNFα and IL1β in monocytes, thereby promoting the inflammatory cascade.

Some metabolites associated with amino acid metabolism differed significantly between the high- and low-PAH exposure groups. Similar alterations were observed in BaP-exposed lung cancer H460 cells^[31]. Hippurate is a product of the conjugation of benzoate with glycine, which leads to an increase in hippuric acid levels and a subsequent reduction in glycine levels. This decrease in glycine levels results in decreased urea and uric acid excretion in the urine^[38]. Consistent with our decrease in urea, a decline in uric acid concentration has also been reported among children around 7 years of age and adults residing near coking plants with high environmental PAH concentrations^[30]. Given that our study collected samples from individuals aged between two and three years old, excessive intake of benzoate preservatives through beverages can be disregarded. Benzoate primarily originates from the microbial degradation of dietary aromatic compounds, such as polyphenols, purines, AAAs, and other amino acids^[39]. In other words, variations in urinary hippuric acid levels are predominantly influenced by benzoate generation via intestinal microbial metabolism of dietary aromatic compounds^[40]. Consistent with the changes observed in hippuric acid levels, gentisic acid also showed a positive association with PAHs. Gentisic acid is a metabolite of the tyrosine metabolic pathway that increases concomitantly with an elevated Lactobacillus abundance^[41]. Additionally, studies have reported that increased Clostridia spp. along with decreased Bifidobacterium levels can lead to higher hippuric acid levels^[42]. Clostridia spp., Bifidobacterium, and Lactobacillus are involved in tyrosine metabolism^[41]. Therefore, PAHs in the human body can affect the metabolism of amino acids, including AAAs. Alterations were also detected in adenosine, a cAMP pathway metabolite. Activation of the cAMP pathway occurs when there is disruption of the glutamate receptor (NMDAR)^[43]. Transcriptomic analysis has revealed that PAH-induced neurotoxicity is more likely caused by interference with neurotransmitter receptors, particularly NMDAR subunits, rather than by genotoxicity^[44]. As neuroactive ligands, AAAs typically enter nerve terminals through receptor-mediated mechanisms to facilitate neuronal growth, development, and functional integrity^[45]. Differential metabolites such as hippuric acid and gentisic acid are involved in phenylalanine and tyrosine metabolism, respectively, both of which are AAAs. Reduced levels of these neuroactive ligands can hinder the production of various types of receptors, such as NMDAR, leading to the disinhibition of glutamatergic output from neurons that inhibit GABA signaling. Consequently, this disrupts the central signal-to-noise ratio, ultimately resulting in abnormal neurobehavioral manifestations and cognitive impairment.

Glycerophosphocholine concentration was positively associated with naphthalene and phenanthrene metabolites. Both naphthalene^[46] and phenanthrene^[47] inhibit acetylcholinesterase activity, leading to accumulation of the neurotransmitter acetylcholine (ACh) and an increase in glycerophosphocholine levels. Moreover, occupational exposure to PAHs is linked to cognitive impairment and elevated serum ACh levels^[48]. ACh is present in both the central and peripheral nervous systems and plays a role in cognitive function, particularly in memory and learning. Interestingly, we also observed that the pantothenate levels were significantly higher in the high PAH exposure group than in the low PAH exposure group. Pantothenate plays an essential role in Ach synthesis^[49]. Exposure to PAHs induces neurodevelopmental toxicity by influencing ACh levels, thereby disrupting glycerophosphocholine metabolism.

Exposure to environmental chemicals induces inflammation in the body, leading to alterations in the intestinal flora and subsequent increases in the concentrations of choline and its compounds, such as glycerophosphocholine^[50]. This suggests that PAH exposure may elicit a systemic inflammatory response, disrupt the composition of the intestinal flora, and further elevate choline levels within the body. The concentrations of deoxyinosine, hypoxanthine, N4-acetylcytidine, and 2-methylguanosine were significantly positively correlated with the abundances of Subdoligranulum and Coprococcus. The genera Subdoligranulum and Coprococcus, members of the Lachnospiruceae family or Clostridial clusters, are known for their anti-inflammatory effects and contribution to intestinal health through components such as butyrate^[51]. They are also associated with neurodevelopmental and neuropsychiatric disorders^[52]. In addition to an increase in deoxyinosine, hypoxanthine, and xanthine levels, we observed a decrease in urea levels. Urea is a downstream metabolite of uric acid. Notably, one-third of uric acid in humans originates from the intestinal flora, such as Escherichia coli of Proteobacteria, which can secrete xanthine oxidase^[53], an enzyme that converts hypoxanthine to xanthine and further metabolizes it into uric acid. Interestingly, we observed a positive correlation between the Proteobacteria genus Haemophilus and increased levels of pseudourine and N4-acetylcytidine, two nucleosides associated with tRNA post-transcriptional modifications. Haemophilus influenza is possesses pseudouridine synthase enzyme that is responsible for converting uridine into pseudouridine at specific positions in the ribosomal RNA^[54]. Among the gut microbiota, the family Comamonadaceae, genus Thalassospira, and genus Burkholderia of the phylum Proteobacteria are positively linked to autism spectrum disorder (ASD), whereas the genus Haemophilus is more prevalent in oral samples from children with ASD^[52]. Pantothenate is an essential nutrient for the growth of certain pantothenate-non-producing bacteria, such as Firmicutes, within the intestine and can enhance ACh synthesis during Lactobacillus proliferation^[49]. Genomic analysis has revealed limited de novo synthesis of pantothenic acid in Bacteroidetes and Proteobacteria^[55].

The altered metabolite, L-gulonic gamma lactone, was positively correlated with Fusobacterium and Flavonifractor. Fusobacterium has been linked to various inflammatory conditions such as bowel disease, atherosclerosis, organ abscesses, adverse pregnancy outcomes, and colorectal cancer^[56]. Similar to Subdoligranulum, Flavonifractor produces acetic and butyric acids, while promoting intestinal health by reducing inflammation^[57].

PAHs promote acetate accumulation, primarily derived from substrates hydrolyzed by acetyl-CoA, which is converted to acetyl phosphate and then to acetate^[58]. Our findings suggest a positive association between Coriobacteriales from the Actinobacteria phylum and acetyl phosphate. The Actinobacteria phylum plays a crucial role in maintaining gut homeostasis, with its involvement in modulating gut permeability, immune system function, metabolism, and the gut-brain axis^[59]. Additionally, Bifidobacteria, a widely used probiotic within the Actinobacteria phylum, has shown beneficial effects in various neurodevelopmental disorders^[59].

The major limitation of this study was its small sample size; therefore, caution should be exercised when interpreting the results. Owing to the limited sample size, it may not have been feasible to detect significant differences in child development and neurobehavioral tests between the high- and low-exposure groups. The small sample size might have resulted in low statistical power and hindered further multivariate analysis, exploration of potential confounding factors, and the exclusion of possible chance findings. Hence, we refrained from performing a multivariate analysis and presented the correlation coefficients using a heatmap and hierarchical clustering technique. The heatmap provides an intuitive visualization of the data distribution, with different colors representing high, low, and middle values for the quick identification of patterns and trends that aid in identifying correlations within the dataset. Second, another notable limitation is our failure to consider potential confounding chemicals, such as phthalate esters, polybrominated diphenyl ether, and polychlorinated biphenyls found in settled house dust, which may affect neurobehavior in children. Future studies should consider the possible combined effects. Third, we only measured average PAH concentrations over a 5-day period at approximately three years of age, potentially overlooking seasonal variations in exposure levels. Finally, the differential metabolites identified using our non-targeted metabolomics approach were relatively quantitative, primarily serving as a means to discover potential differential metabolites and generate hypotheses^[21]. These differential metabolites require further quantitative verification in pediatric populations.

Heightened exposure to PAHs in young children may disrupt their overall metabolic equilibrium, particularly affecting specific metabolites originating from the gut microbiota. This can lead to variations in metabolites, including amino acids, lipids, carbohydrates, nucleotides, and vitamins, across various metabolic pathways. Although these differential metabolites could distinguish group differences, their effectiveness in predicting individual exposure was limited. Moreover, a significant proportion of these differential metabolites and associated metabolic pathways were linked to inflammatory responses and were potentially involved in modulating the intestinal flora, thus potentially influencing neurodevelopment. Owing to our small sample size limitation, it is crucial to increase the sample size and conduct further research to quantify PAH exposure for clinical and public health applications. Additionally, gaining a deeper understanding of the molecular mechanisms underlying the neurodevelopmental abnormalities caused by PAHs is essential, with a specific focus on the involvement of the gut-brain axis and the regulation of intestinal flora.