The HERACLES Greek cohort focused on the impacts of exposure to heavy metals (originating from a waste disposal site) on child neurodevelopment in Athens, Greece. The study was initiated in 2012, and 300 children aged 3 to 8 were enrolled, who lived in the proximity of the waste management site (landfill, children with a residential address up to a distance of 12 km were enlisted). Characteristics of the study population including demographics and health outcomes are given in Table 2.

For UPLC-HRMS analysis, urine and plasma samples were thawed under stable conditions to standard protocols described from Want et al. (35) and Theodoridis et al. (36).

Urine: 600 μl from urine samples was centrifuged at 10,000 rpm for 10 min. A supernatant of 500 μl was placed on autosampler vials and diluted with 1:2 of LC-MS water. After that, the samples were ready for analysis. The autosampler operated at 4°C.

Plasma: A quantity of 200 μl plasma sample was transferred to a new Eppendorf and diluted with 600 μl of cold methanol. After that, the samples were centrifuged at 10,000 rpm for 20 min. A supernatant of 300 μl was placed to a clean Eppendorf and dried using a Techne Sample Concentrator. All samples were reconstructed with 100 μl of LC-MS water. The samples were further centrifuged (10,000 rpm for 20 min). A final quantity of 95 μl of supernatant was placed into 2-ml vials with inserts and placed in UPLC autosampler operating at 4°C, ready for analysis. In addition, a 50 μl quantity was used for the pooled quality control (QC) sample.

Sample analysis was performed on a ThermoFisher Scientific model LTQ Orbitrap Discovery MS, with a resolution of 30,000. The spectra from both urine and plasma samples were acquired in both positive and negative ionization modes. The mass scanning range was set at 50–1,000 m/z. The capillary temperature was set at 320°C. Nitrogen sheath gas and auxiliary gas was flow rate was set to 40 L/min and 8 L/min respectively and the spray voltage at 4.5 kV. LC-MS uses a gradient of two solvents. One hundred percentage LC-MS water with 0.1% formic acid as mobile phase A and 100% methanol with 0.1% formic acid as mobile phase B. For urine samples analysis, the flow rate was 500 μl/min for both the positive and negative modes. For the plasma sample analysis, the flow rate was 400 μl/min in positive mode and 360 μl/min in the negative mode.

For chromatographic separation, an Acquity UPLC HSS T3 column (100 × 2.1 mm, 1.8 μm, Waters, Milford, MA, USA) was used, which was kept at a constant temperature of 40°C. The gradient of the mobile phase for urine samples for both positive and negative mode was: 1% B at 0 min, 1% B at 1 min, 15% B at 3 min, 50% B at 6 min, 95% B at 9 min, 95% B at 10 min, 1% B at 10.1 min and 1% B at 14 min. For the plasma samples, a different gradient was used; for the positive mode: 0% B at 0 min, 0% B at 1 min, 100% B at 16 min, 100% B at 20 min, 0% B at 22 min, 0% B at 24 min. For the negative mode: 0% B at 0 min, 0% B at 2 min, 100% B at 17 min, 100% B at 22 min, 0% B at 24 min, 0% B at 26 min. Blank samples, and pooled QC samples, were used to monitor the systematic signal deviations between the batches. Two blank samples containing internal standards in known concentrations were injected before each batch, in the middle and at the end of each batch, for (a) checking the condition of the column and (b) calculating the mass error in ppm of the isotope patterns. Additional, ten QC pooled samples were injected at the onset of the experiment for the conditioning of the column of the LC system. The pooled QC sample was repeated every ten samples in urine samples analysis, whereas plasma samples the QC samples were repeated every five plasma samples.

First samples were centrifuged at 10,000 rpm for 10 min. Then, 500 μL of the supernatant was transferred to a new Eppendorf and mixed with 120 μL of buffer solution (Na₂HPO₄,0.2 M, NaH₂PO₄,0.3 Min 50% D₂O/50% H₂O) and 0.1% TSP-d4, which was used as chemical shift reference (δH 0.00 ppm). The pH of the samples was 7.337. All samples were vortexed and placed in −4°C for 7 min. After 7 min, samples were thawed and centrifuged, and the supernatant was transferred to a 5 mm NMR tube. A 1 ml aliquot was evaporated to dryness under vacuum for plasma samples. Reconstitution was performed with 50 μl deuterium oxide (D₂O), followed by vortex and centrifugation at 14,000 rpm for 5 min at 4°C. The supernatant was evaporated once more and then reconstructed in 660 μl of 100 mM phosphate buffer containing 1 mM trimethylsilyl propanoic acid (TSP) and 1 mM sodium azide. The extract was vortexed before a final centrifugation step at 14,000 rpm for 5 min at 4°C. The supernatant was transferred to a 5 mm NMR tube.

For urine samples analysis, a 600 MHz Varian spectrometer was used. The spectrometer frequency was 599.938 MHz in an OneNMR Probe and a ProTune System (Agilent) using on-resonance pre-saturation to suppress the intensity of the water signal. Proton chemical shifts typically range from −2 to 10 ppm (spectral width 9615.4 Hz), with 128 scans, a relaxation delay of 2 s, acquisition time 4 s and pulse width 8.587 μs. Plasma samples were acquired in a Bruker Avance 500 MHz NMR spectrometer equipped with a TCI cryoprobe. Data acquisition and processing were performed using the software package Topspin v 1.3 (Bruker, Germany). The central frequency used was 500.1323505 MHz, using on-resonance pre-saturation to suppress the intensity of the water signal, followed by a 1D NOESY pulse sequence with irradiation of the residual water signal the mixing time (200 ms). The observation pulse length was set at 10.0 μs, the delay between transients was 3 s, and 65,536 complex data points were acquired with a spectral width of 10,400 Hz (corresponding to a chemical shift range of 14.0019 parts per million, ppm), giving a final acquisition time of 4.679 s. The total time of the experiment was ~67 min and included eight unrecorded (dummy) transients and 512 acquisition transients (scans).

Raw data generated from positive and negative ionization were pre-processed as two different experiments. The tool msConvert included in the ProteoWizard toolkit (37, 38) was used to translate the data into the.mzML open format. Spectral processing was performed using the Bioconductor R—based packages XCMS v.3.10.1 (39) running under R version 3.6.1 (https://www.r-project.org/). Chromatographic peak detection was performed using the centWave algorithm, and the Obiwarp method was used for alignment. The definitions function was used for peak correspondence. We used the fillChromPeaks method to fill in intensity data for missing values from the original files due to false negatives. Finally, the PerformPeakAnnotation function is used for isotope and adduct annotation using the CAMERA package (40). The resulting matrix was further reduced by the 80% rule applied to the QC samples to obtain consistent variables. The instrument and overall process variability were then determined by calculating the median RSD for authentic internal standards and all endogenous metabolites. Normalization by the median, mean centering scaling, and log transformation was performed to transform the data matrix into a more Gaussian-type distribution, thus reducing systematic error in experimental conditions. The annotation R package xMSAnnotator (41) was used to perform the accurate mass carries in online compound databases [HMDB (42), LipidMaps (43), and KEGG (44)]. The adduct list used for database matching included “M+2H,” “M+H+NH4,” “M+ACN+2H,” “M+2ACN+2H,” “M+H,” “M+NH4,” “M+Na,” “M+ACN+H” in positive ionization mode, and the following in negative ionization mode: “M+ACN+Na,” “M+2ACN+H,” “2M+H,” “2M+Na,” “2M+ACN+H” “M-2H,” “M-3H,” “M-H2O-H,” “M+Na-2H,” “M-H,” “M+Cl,” “M+FA-H,” “M+K-2H,” “2M-H.” The list of the detected candidate metabolites was filtered by metabolite status (Detected) and biospecimen (Blood). Confirmation of identified biomarkers was performed by comparing the RT and fragmentation pattern of authentic analytical standards from the in-house library or MS/MS spectra available in databases like HMDB and Metlin (45).

Spectral analysis of NMR data proceeded using MestReNova (Mnova 11.0.3) (http://mestrelab.com), while for metabolite identification ChenoMx (http://www.chenomx.com) was used in addition. Briefly, after loading the spectra of all the samples using the superimposed command, the first step was to correct the position of the reference peak sample. The detected peaks (binning) were grouped using length values lower than 0.04 ppm in this study; the used reference was deuterium oxide (D₂O) due to the used buffer. After the reference correction, the alignment of all reference peaks was checked. The Smoother Whittaker algorithm was chosen for baseline correction; after that, the spectrum phase was checked and corrected. For phase correction, an automatic algorithm is preferable. The detected peaks (binning) were grouped using length values lower than 0.04 ppm. Then, the spectrum was imported into the ChenomX NMR Suite 8.2 for peak identification. Finally, the previously identified peaks of TSP and metabolites were integrated using the MNOVA software, and the metabolite peaks were identified. In cases where multiple peaks and peaks characterized a metabolite in different areas of ppm, the area of all these peaks were added to fill the corresponding cell on the sheet of import file to MPP.

All LC-MS and NMR data were deposited to the EMBL—EBI MetaboLights database with the MTBLS1882 identifier.

Pathway analysis has been carried out by combining two GeneSpring modules, Mass Profiler Professional (MPP) and Pathway Architect (Agilent Technologies). To calculate the probability of enriching the listed metabolites with a pathway, a hypergeometric test was used, while the p-values have been corrected for multiple testing using the Benjamin-Hochberg method (46).

The exposome-wide association study (EWAS) paradigm has been coined Patel et al. (11) to describe the associations among several variables composing the human exposome, using unsupervised learning (47). The study predictors can be classified in numerical, nominal categorical, or ordinal categorical. For example, the exposure factors were continuous variables since they were data from mass spectrometry. In contrast, the variables regarding the frequency of specific food items or water consumption were ordinal. Before applying the logistic regression algorithm, we investigated the skewness of the continuous variables and then transformed them into logarithmic form. Next, we adjusted each observation to the mean and scaled it by the standard deviation and after that it was classified into bins. The number of selected bins is 10. In the present study the range of each class is calculated taking the difference between the maximum and minimum value divided by the number of bins. Therefore, observations are clustered /into categorical classes so that they can then be used in logistic regression. The rest of health outcome variables that were used are presented in the following table in order to simply describe the limits of each bin. In addition, we converted categorical variables into numerical ones, using integer encoding, whereby each unique label was mapped to an integer (Table 3).

Henceforth, we used survey-weighted logistic regression models to associate each of the exposure factors and modifiers with a health outcome like the psychomotor scores (cognitive, language, and motor score), while adjusting for child gender, and a variety of sociodemographic, anthropometric and other gestational and post-delivery parameters collected from self-reported questionnaires. Spearman correlation, a non-parametric test, was applied to calculate the correlations among variables and avoid any distributional assumptions. The False Discovery Rate (FDR) q-value was calculated to associate the factor with a health outcome, like Bayley test levels, controlling the type I error using the Benjamini-Hochberg step-down approach (48). Applying permutation resampling to Bayley III test scores allowed us to validate the FDR results (49). The significance level was set at 0.05, corresponding to an FDR of 10%. The same procedure was repeated for the study of exposure factors and the detected metabolic pathway associations, as well as for revealing the associations between health outcomes (Bayley III test) and the perturbed metabolic pathways. An Ubuntu 16.04 server was used for carrying out the calculations, and the “X-Wide Association Analyses” R package (50) was used for the logistic regression and the FDR calculations, while the “RCircos” package (51) was employed for visualizing the results through correlation globes. Dataset skewness correctness has been checked using the R library “moments,” while the R package “permute” was used for the restricted permutations of the randomization tests. Finally, for data analysis, clustering, utility operations, computing sample size and power, the Hmisc R package (52) was used to import and annotate datasets, and impute missing values.

In order to generate a list of candidate biomarkers for the stressors of interest, which were identified by untargeted metabolomics analysis, a manual literature search was carried out including combinations of the following terms: “exposure,” “metals,” “neurodevelopment disorders,” “biomarkers” and “underlying mechanisms.” The query was performed for both levels of analysis (molecular biomarkers and pathway analysis) for all the experiments. It was designed not only to determine candidate biomarkers based on their known biological function but also to avoid the exclusion of biomarkers associated with the occurrence of a stressor, even though there is yet lack of knowledge regarding its exact biological function.