Introduction

Psychological stress has systemic effects on the body that may lead to chronic impairments in physical and mental health, resulting in reductions in longevity and quality of life1,2. More than 70% of human diseases are considered to be stress-related3, with the most common stress-related disorders being posttraumatic stress disorder (PTSD), anxiety disorders and major depressive disorder (MDD)4. Despite the established importance of stress in disease, we currently lack objective measures of the biological impact of stress on the body. A better understanding of the metabolic changes caused by stress is crucial for elucidating the mechanisms underlying the development of stress-related diseases and developing better treatments3. Furthermore, personalised medicine requires biomarkers that are sufficiently sensitive and specific for diagnosing potentially damaging stress exposure and predicting future outcomes5.

Elevated cortisol, which can be measured non-invasively in saliva or hair, is one of the most common measures currently used to assess acute and chronic stress. However, due to the inverted U-shaped relationship between cortisol levels and numerous neurobiological and behavioural endpoints, cortisol is nonspecific and cannot be used to reliably distinguish understimulation from overstimulation and stress6. Moreover, there are individual differences in the level of stress associated with optimum stimulation for health and well-being with individuals vulnerable to stress-related disorders having lower optimum cortisol levels than individuals resilient to stress6.

Due to the complexity of the stress response, a more holistic approach is required, moving away from the measurement of single biomarkers such as cortisol. Triangulation, whereby multiple biomarkers, each with its own strengths and weaknesses, are measured simultaneously, should deliver a more reliable estimate of an individual’s condition. Untargeted metabolomics identifies and compares the relative concentrations of a large cohort of metabolites (small molecules < 1500 Da) between samples of interest, yielding a specific metabolic signature or fingerprint of a given exposure or manipulation7,8,9. The results can be used to identify a panel of metabolites that effectively discriminate against different conditions. We test the hypothesis that untargeted metabolomics conducted on samples obtained before and after exposure to a major stressor will reveal new biomarkers and shed light on the metabolic pathways involved in the etiology of stress-related diseases.

Collecting samples longitudinally before and after a terrifying, life-threatening event, characteristic of a major stressor, is practically and ethically challenging. Here we make use of an opportunistic pig model that we argue has strong ecological validity for the types of events known to potentially result in PTSD and other stress-related disorders. Handling of pigs pre-slaughter involves a 24 h period during which animals are transported, re-mixed with unfamiliar pen-mates and held in a new environment where they are exposed to the noise and smell of the slaughterhouse10. These events have been widely described in the literature as extreme stressful events for pigs10,11,12,13. These events result in prolonged thirst and hunger, restricted movement, negative social behavior, resting problems, fatigue, pain, fear and respiratory distress12, with physiological consequences sufficient to impair meat quality10,13. Pigs have physiological similarities to humans making them a common animal model in biomedical research14,15. Additional benefits of the model include the homogeneous genetic background, age, nutrition and rearing environment of the animals, all of which reduce random variation in physiology. Saliva samples were collected non-invasively from standard environmental enrichments provided in the animals’ pens, ensuring that baseline samples were uncontaminated by sampling stress.

Results

Effects of 24 h of stress on the saliva metabolome

For the examination of acute stress as detailed in Fig. 1, longitudinal saliva samples were collected non-invasively from 200 pigs by providing cotton ropes for the animals to chew as environmental enrichment16,17,18. The pigs were penned in groups and each rope provided a pooled sample from 5 to 10 animals. Saliva samples were collected from the pigs twice: (1) 24 h before slaughter, while pigs were still in their pen, in a familiar group from weaning (n = 31 samples); and (2) after short transport, regrouping with unfamiliar pen-mates, and 24 h holding in the new environment of the slaughterhouse (n = 32 samples). All samples were kept on ice, and stored at -80 °C within two hours of collection until the metabolomics analyses. An untargeted metabolomics approach was used involving ultra-high-performance liquid chromatography coupled with high-resolution mass spectrometry (UHPLC-HRMS). Metabolites were extracted from the saliva samples using an extraction procedure (methanol:acetonitrile) suitable for polar metabolites. UHPLC-HRMS was performed using a hydrophilic interaction liquid chromatography (HILIC) analytical column, in both positive and negative ionization mode for broader metabolite coverage. A pooled sample from all samples was injected from the same vial for quality control (QC), as well as pooled samples extracted four times for quality assurance (QA). Overall, obtained UHPLC-HRMS data revealed the relative quantification of 2518 metabolite features in total, in the positive ionization mode, and 1165 metabolite features in the negative ionization mode. After data evaluation and cleaning, including background spectral filtering (noise elimination) and QC based filtering (metabolite features with coefficient of variation (CV) ≥ 30% were removed), 1564 metabolite features from the positive mode, and 850 metabolite features from the negative mode remained for further statistical analyses. Of these, 920 metabolite features were putatively annotated in the positive ionization mode, and 531 in the negative ionization mode.

Figure 1
figure 1

Study design. For the examination of acute stress, saliva samples were obtained non-invasively from 200 pigs. Cotton ropes were provided, yielding 63 pooled saliva samples, each from 5 to 10 pigs. Saliva was collected from the same pigs twice: (1) before stressful conditions (n = 31); and (2) after 24 h in stressful conditions (n = 32). Metabolites were extracted and analysed using an ultra-high‐performance liquid chromatography high-resolution mass spectrometer. Raw data was used for data pre-processing, identification of features as metabolites, and data analyses for the identification of altered metabolites and metabolic pathways. Procreate (version 5.2.6) and Adobe Photoshop Creative Cloud (version 23.5.1) were used for the graphics and illustrations.

Full size image

A principal component analysis (PCA) score plot revealed clear separation between pig saliva metabolites before and after the acute stress intervention, suggesting a clear influence of the intervention on the saliva metabolome (Fig. 2).

Figure 2
figure 2

PCA score plots of pig saliva metabolome. A Principal Component Analysis (PCA) plot of the metabolome data that characterizes pig saliva before (Blue) and after (Orange) acute stress. Each dot represents a sample and the colour represents the type of sample. Repeatability (Grey) is demonstrated by pooled samples, injected along the analyses as quality controls. (a) A PCA score plot of metabolite features from data acquisition in positive ionization mode. (b) A PCA score plot of metabolite features from data acquisition in negative ionization mode. Graphs were created by GraphPad Prism (version 9.4.).

Full size image

Furthermore, relative concentrations of 932 metabolite features in the positive ionization mode, and 591 metabolite features in the negative ionization mode were significantly different before and after the stressful intervention (Fig. 3; Adj P < 0.05). Of these features, those with at least a two-fold change in relative concentration were considered differentially affected by stress. Accordingly, in the positive ionization mode 243 metabolite features were significantly upregulated (a fold change of 2.00 to 17.12) after stress, and 357 metabolite features were significantly downregulated (a fold change of 2.00 to 333.33) after stress (Fig. 3a; Adj P < 0.05). In the negative ionization mode, 49 metabolite features were significantly upregulated (a fold change of 2.01 to 53.05) after stress, and 333 metabolite features were significantly downregulated (a fold change of 2.04 to 250.00) under stress (Fig. 3b; Adj P < 0.05).

Figure 3
figure 3

Differential metabolite features after stress. (a–b) Differential analysis of metabolite ratios between pigs in their familiar environment (Before) and 24 h after (After) the onset of stress. (a) Volcano plot of metabolite features from data acquisition in positive ionization mode. (b) Volcano plot of metabolite features from data acquisition in negative ionization mode. Each dot represents a metabolite feature. Significant upregulated metabolite features are in Red on the right (2 \(<\) fold change; Adj P < 0.05). Significant downregulated metabolite features are in Blue on the left (fold change \(<\) 0.5; Adj P < 0.05). (c–d) Heat map analysis of differential metabolite features (fold change \(<\) 0.5; 2 \(<\) fold change; Adj P < 0.05), (c) in both positive and (d) negative ionization modes. Each row represents a sample and each column a metabolite feature. Purple indicates that metabolite features were downregulated, and Yellow indicates that metabolite features were upregulated. Graphs a, b were created by GraphPad Prism (version 9.4.); Heat maps c, d were plotted by R (version 4.2.0), using R package pheatmap (version 1.0.12).

Full size image

Follow-up analyses included investigation of metabolite features responsible for the clustering of the metabolome changes before and after the stressful intervention. The complete screening for biomarkers is summarized in Table 1, and fully detailed in Table S1. As mentioned, we searched for metabolite features for which the relative concentrations after stress were not only significantly different, but also showed at least a two-fold change (Adj P < 0.05; 2 ≤ FC or FC ≤ 0.5). Each of these metabolite features that could serve as potential biomarkers was manually checked for the match of the precursor ion mass (m/z) and fragments (MS/MS) to online and in-house libraries. Metabolite identification confidence level for each identified metabolite in this study was assigned following the modified four-level classification scheme from the Metabolomics Standards Initiative19,20: level 1- match to commercial standards; level 2- full match to online databases of the precursor ion and its fragments (both m/z and MS/MS); level 3- putative annotation based on precursor ion; and level 4- not annotated). To ensure higher confidence, the metabolites discussed as potential biomarkers in this paper are only those with identification confidence levels 1 and 2. These were mainly decreased relative abundances of amino acids and metabolites synthesized endogenously by amino acids and B vitamins. In addition, the most promising biomarkers that were not identified with high confidence, or were not annotated at all, are detailed in Table S1.

Table 1 Summary of potential biomarkers.
Full size table

Another approach in the investigation of novel biomarkers and metabolic fingerprints, is to look at the alteration of entire metabolic pathways, rather than single potential biomarkers. Therefore, 144 metabolite features from both the positive and negative ionization with higher annotation confidence (level 1 and level 2; when both precursor ion and its fragments, MS/MS, were matched to in-house library of commercial standards or online databases). 95 of these putatively annotated metabolite features were successfully matched to a metabolic pathway by Metaboanalyst software 5.021. Several metabolic pathway alterations were induced following the 24 h stressful intervention. These pathways are represented by both −log10 (Adj P-value) \(>\) 1.1220, which is equivalent to a Adj P-value < 0.05 and high pathway impact score, calculated by overrepresentation of metabolites from the same pathway, as well as highly important metabolites to the pathway. Accordingly, as shown in Fig. 4, histidine metabolism; phenylalanine, tyrosine and tryptophan biosynthesis; arginine and proline metabolism; and arginine biosynthesis are significantly different, as well as have relatively high pathway impact. Moreover, the pathways of nicotinate and nicotinamide metabolism; lysine degradation; alanine aspartate and glutamate metabolism; and vitamin B6 metabolism were also highlighted and show a promising direction for further investigation.

Figure 4
figure 4

Pathway analysis plot of altered metabolic pathways after acute stress. Dots represent metabolic pathways altered after 24 h of major stress. The Y-axis represents the logarithm transformed P-value adjusted for multiple comparisons. The X-axis represents the pathway impact, calculated by overrepresentation of metabolites from the same pathway, as well as highly important metabolites to the pathway. The pathways that were most associated with stress are a combination of high pathway impact and -log10 (Adj P-value) \(>\) 1.1220, which is equivalent to a P-value < 0.05 (the top right corner of the graph). Colour represents significance, from White (not significant) to Red (significant), and larger dot size represents a higher impact score. The figure was created by the online tool MetaboAnalyst 5.0 (www.metaboanalyst.ca/MetaboAnalyst/).

Full size image

Discussion

The trajectory of stress and related disorders has been widely investigated1,2,3, but there remains a gap in our understanding of the acute systemic effects of major stress3,22,23. Our aim was to characterize the metabolic fingerprint of major stress, in order to provide a comprehensive understanding of its acute physiological effects and to reveal novel potential biomarkers. We utilized an untargeted metabolomics approach in which we used UHPLC-HRMS to characterize the change in the saliva metabolome of the same group of pigs measured before and 24 h after exposure to major stress involving transport, social mixing and holding in the novel and frightening environment of a slaughterhouse. Our results provide a basis for further investigation of the impacts of stress and the mechanisms underlying the development of stress-related disorders such as PTSD in humans. Moreover, they may shed light on potential drugs and nutritional supplements for treating stress-related disorders.

In this study, we found several hundreds of metabolite features that display a 2.00 to 333.33 fold change after major stress (600 metabolite features in the positive ionization mode and 382 metabolite features in the negative ionization mode). These metabolite features (detailed in Table 1 and Table S1) are divided into two main groups; The first, are metabolite features for which the retention time, as measured by the liquid chromatography, and the molecular mass of the molecule and its fragments, as measured by the mass spectrometer, could be used to identify the feature as a specific metabolite (identification confidence levels 1 and 2). The second group, are metabolite features that are part of the metabolic fingerprint of stress and may provide potential biomarkers, but for which identification requires additional analytical methods (confidence levels 3 and 4). The main metabolic changes identified were decreases in amino acids and metabolites synthesized endogenously by amino acids and B vitamins. Many types of stress including psychosocial stress, injuries, surgeries, infection and cancer activate metabolic pathways that consume amino acids24,25,26. In this study, the ratio of relative concentrations of arginine after versus before stress was 0.154, indicating that arginine was significantly lower after 24 h of stress. Arginine is a conditionally essential amino acid, because while it can be synthesized, it becomes essential during trauma and disease24,27,28. Lysine and tryptophan also decreased after stress in this study (ratios of 0.162 and 0.492 respectively). Prolonged lysine insufficient diets are correlated with impaired quality of life and mental health29. Thus, a decrease under prolonged stress might be consequential, especially if the diet is deficient in lysine. Tryptophan is also an essential amino acid, which functions as a precursor of serotonin, as well as in the synthesis of kynurenine. The tryptophan-kynurenine pathway is associated with cognitive symptoms and neurodegenerative diseases30. While the role of serotonin in regulating mood and cognition is known, the complete mechanism of kynurenine and the unwanted symptoms of neurodegenerative diseases is not fully discovered30. There is controversial evidence in the literature regarding the influence of tryptophan and its supplementation in diet on mood change and stress-related disorders such as major depression31,32. B vitamins relative abundances also decreased in our study; the ratio of relative concentrations of nicotinamide after stress versus before was 0.498. Nicotinamide, the amide form of vitamin B3, has a role as a precursor in the synthesis of coenzymes nicotinamide adenine dinucleotide (NAD+) and nicotinamide adenine dinucleotide phosphate (NADPH), which are responsible for various functions including oxidative deamination and lipid catabolism33. Nicotinamide has a key role in development, growth and maintenance of the central nervous system and is also considered to have protective effects on the nervous system, providing protection from neurodegenerative diseases such as Alzheimer’s, Parkinson’s and age-related macular degeneration (AMD)34. Another B vitamin that was relatively decreased in this study after stress was pyridoxine, vitamin B6 (0.359). Pyridoxine has a key role in amino acids metabolism and cell functioning35. Moreover, as for other B vitamins, pyridoxine is considered to have protective effects against cognitive decline and the related pathologies such as Alzheimer’s disease36. The decrease in these three metabolites (tryptophan, nicotinamide and pyridoxine) under acute stress may be important in elucidating the effect of stress on the later trajectories. Additional metabolites that may shed the light on these trajectories, were also decreased under stress in this study. Carnosine decreased after stress (ratio of 0.467). Carnosine, a dipeptide, which is endogenously synthesized from the amino acids alanine and histidine, has many physiological functions and is dominant in tissues such as the brain, heart, muscles, liver, brain, and kidneys37. L-carnosine has been shown to function as an antioxidant, an anti-inflammatory, an accelerator of wound healing and it might even have antidepressant effects38. Furthermore, L-carnosine has been shown to act as an inhibitor of the toxic amyloid-β peptide (Aβ) present in the vascular plaques in Alzheimer's patients, as well as in reducing aging-related mitochondrial dysfunctions37,38,39,40. Hypoxanthine which also decreased under stress (0.416) is a derivative of purines, and is associated with abnormal purine metabolism in Alzheimer’s patients41. Agmatine, which is also induced endogenously as part of the stress and inflammatory response, has a role as a neuromodulator. It has been suggested that agmatine has a protective effect on the development of trauma and stress related disorders, such as PTSD, anxiety, and depression42. Moreover, it may have protective effect against neurodegenerative diseases such as Alzheimer's and Parkinson43,44,45. In this study after 24 h in stressful conditions, saliva agmatine also decreased (ratio of 0.377). It has been shown that as a result of stress, trauma and inflammation, agmatine level is increased, but cannot compensate fully for its harmful impact, at least not the endogenous production46,47. We hypothesize that the agmatine that was produced under stress was consumed as part of the protective mechanism for neural functioning, and therefore was already decreased after 24 h in stressful conditions. It has been described that agmatine is elevated after the first few hours from stressful exposure, but the exact mechanism of agmatine as part of the stress response, as well as the longer term effect remained unclear48. Moreover, it was also demonstrated in rodents that agmatine has anxiolytic effect, as well as suppression of depression-like behaviour49,50. Further investigation requires measuring the levels in saliva and other biofluids, as well as in the brain, from the stressful exposure, and after different time periods, in order to get a comprehensive understanding. Individual sampling might also be considered in complementary future research. In this research, group level sampling was chosen in order to ensure results are not affected by the sampling method. Individual sampling is stressful for the animals and could not be performed non-invasively, without direct contact with the pigs in the study. This is a trade-off between group sampling to non-invasive sampling. Here we decided it is best to compromise on group samples rather than compromising on the sampling method. Thus, we ensured that baseline samples were uncontaminated by sampling stress.

The main limitation of this study, as in all untargeted metabolomics studies, is the inability to identify several metabolite features. Identification of metabolites is the main bottleneck in the workflow in untargeted metabolomics. Identification is based on the mass of the intact metabolite (more specifically, the mass to charge ratio, m/z), its fragments (MS/MS), and the retention time in the analytical column. Although we can record this information for thousands of metabolite features in one analysis, still many of the metabolite features remain unidentified as a specific metabolite. This limitation is due to the lack of available standards, limited databases and the available algorithms51. It is likely that several unidentified metabolites have the potential to serve as novel biomarkers and are thus worthy for further investigation. For example, Table S1 includes metabolite features that are significantly differentially regulated after stress, with at least a ten-fold change ratio (Adj P < 0.05; FC > 10; FC < 0.1). These features are putatively annotated as level 3 or not annotated at all (level 4). One future solution to the large number of unidentified metabolite features will be a cross-species analysis of the effects of stress. Triangulating the overlap in the metabolite features altered by different types of stress, in different species, and in different biological matrices, might help to narrow down the list of potential biomarkers from hundreds, to a shorter list of interest for further identification and validation.

In light of all the above, this research shed light on the metabolic changes under acute stress, by metabolomics analyses in a unique pig saliva model. It lays the foundations for further research for the identification of the unidentified metabolite features, as well as provide a new perspective about stress. Previous studies have focused on upregulated metabolites, such as cortisol, that are elevated under stressful conditions. The current research emphasizes a decrease in several metabolites that are downregulated under acute stress. The decrease in these metabolites, several of which might have protective effects on the nervous system, reducing the likelihood of neurodegenerative diseases and psychological disorders, suggests a novel future direction in research on trauma and stress-related disorders. Perhaps, focusing on downregulated metabolites, as we did in this research, would lead to new future discoveries. This is similar to what has been suggested regarding stress-related disorders, whereby focusing on the mechanism of resiliency, rather than the stress, may bridge this gap of knowledge and lead to new therapeutic targets52,53. Clinically, focusing on these downregulated metabolites as prophylactic treatment may allow individuals to cope better with stress, and stay for longer on the left side of the inverted U-shape, as resilient, before getting to the turning point of vulnerability. There are studies in both humans and animals that provided amino acids and other metabolites that were relatively decreased under stress in our study as experimental supplements. However, probably due to the complexity in these kinds of nutritional studies under specific conditions, as well as the gap of knowledge in the systemic stress response, there is not sufficient scientific-based evidence for the use of these metabolites as treatment. In addition, as far as we are aware, there have been no studies that examined the use of these metabolites prophylactically. Future direction may include targeted analysis of the concentration of these metabolites and other related metabolites for their individual exact concentration before exposure to stress, and prophylactic personal treatment accordingly. Moreover, there are similar opportunistic animal models, including nutritional studies in farm animals that can be beneficial for both human and non-human animals for the research of stress-related diseases and its potential treatments.

In summary, our study indicates that after acute stress of 24 h, in a unique pig model, while several metabolites are upregulated, a higher number of metabolites are downregulated, and their related metabolic pathways changed accordingly. The main metabolic changes that we identified were a decrease in amino acids and metabolites synthesized endogenously by amino acids and B vitamins, as well as potential unidentified metabolite features that require further investigation. The metabolic fingerprint of stress that we reveal provides new insights into the trajectory of stress and the related disorders. This research demonstrates a clear future direction for the investigation of stress and resiliency, as well as potential therapeutics and nutritional targets.

Materials and methods

Sample collection

Saliva samples were collected non-invasively from 200 pigs, ages ~ 165 days, mixed breed of Landrace, Large-White, Pietrain and Duroc. The study was performed at Lahav Animal Research Institute (LAHAV C.R.O; Kibbutz Lahav, Israel) and the non-invasive sample collection protocol was approved by the official Tel Aviv University Institutional Animal Care and Use Committee. All experiments were performed in accordance with relevant guidelines and regulations.

The pigs were housed in groups and each saliva sample represented a pooled sample from 5–10 pigs. Saliva was collected by providing 100% cotton ropes as environmental enrichment, similar to the method recently published16,17,18. In brief, pigs chew the cotton ropes for about 20 min; the saliva is collected by squeezing the rope with a sterile bag, and immediate transferring it into 15 mL plastic tubes; centrifugation of the saliva and aliquots into 2 mL tubes. Saliva was collected from the same 200 pigs twice: (1) 24 h before slaughter, while pigs were still in their pens, in their familiar environment (n = 31); (2) After a short transport to the slaughterhouse (10 min including loading time), random regrouping with unfamiliar pen-mates, and 24 h wait in new environment at the slaughterhouse (n = 32). Although the same 200 pigs are represented in the samples collected both before and after stress, the individual samples represent different combinations of individuals before and after stress due to the regrouping that occurred. Thus, the samples represent independent pooled samples from the same population of pigs collected before and after stress. The samples collected before the stressful conditions are the control for the samples collected after 24 h under stressful conditions. It cannot be ruled out that there are interactions of metabolites with the cotton ropes, but the procedure was exactly the same for samples of the pigs before and after the stressful conditions. The measures are relative abundances and not the concentration in the sample. Here we make use of an opportunistic pig model, taking advantage of the common procedures pigs undergo in the last 24 h of their lives, without any intervention for experimental purposes only. Despite, the ARRIVE guidelines 2.0 were followed54 and the complete checklist was reported. The relevant criteria to this model are reported in the main text.

Saliva samples were stored on ice immediately after collection, and within two hours aliquoted and stored at −80 °C until the extraction procedure and data acquisition.

Sample preparation

Metabolites were extracted from the pig saliva samples using an extraction procedure targeting primarily polar metabolites, similarly to recently published methods55,56. Saliva samples were thawed on ice and briefly vortexed (10 s). Saliva samples were centrifuged at 1300 rpm, 4 °C, for 10 min and the clear fluids were transferred to new 1.5 mL laboratory centrifuge tubes. Samples were vortexed again and 50 μL of each sample, including four QC samples were transferred to a new 1.5 mL laboratory centrifuge tube (from a pooled sampled prepared from 75 µL from each sample in a 15 mL tube as a pooled QC vial. Then, 150 µL of ice-cold 50:50 methanol:acetonitrile (v/v) was added to each. Next, vortex for 15 s and centrifugation at 14,000 rpm, at 4 °C for 20 min was performed. Finally, 150 μL from each sample were transferred to LC vials with inserts for immediate analysis. Samples were injected in randomized order with quality control samples every eight injections.

Quality assurance and quality control

In order to assure the quality of the data, several actions took part: (1) A pooled quality control sample (QC) was prepared by taking 75 µL from each experimental sample, extracted as described above, and injected from the same vial every eight samples. The pooled QC sample injected along the sequence, serves as quality control for the instrument performance. (2) A pooled QC sample was extracted four times (QA), injected randomly along the sequence as quality assurance for the extraction procedure. (3) Extracted water sample was used as a blank, to subtract background ions in data analysis. R package RawHumus57 was used to evaluate the data quality based on QC samples. The resulting QC reports are shown in File S2 for QC data obtained in positive ionization mode and File S3 for QC data obtained in negative ionization mode.

Ultra-high-performance liquid chromatography-tandem mass spectrometry

High performance liquid chromatography was performed by Vanquish UHPLC using analytical Accucore™ 150 Amide HILIC column (100 × 2.1 mm, 2.6 μm), coupled to a high-resolution mass spectrometer Orbitrap Q-Exactive™ Focus (Thermo Scientific™) equipped with a HESI ion source. All samples were first analysed in positive ionization mode and then injected again for analysis in the negative ionization mode. Mass calibration was performed on the day of starting the data acquisition.

Positive ionization mode: mobile phase A: 10 mM ammonium formate, 0.1% formic acid, 95% acetonitrile/water (v/v). Mobile phase B: 10 mM ammonium formate, 0.1% formic acid, 50% acetonitrile/water (v/v). Gradient: t = 0.0 min, 1% B, t = 1.0 min, 1.0% B, t = 3.0 min, 15.0% B, t = 6.0 min, 50.0% B, t = 9.0 min, 95% B, t = 10.0 min, 95% B, t = 10.5 min, 1% B, t = 15.0 min, 1% B. Negative ionization mode: mobile phase A: 10 mM ammonium acetate, 0.1% acetic acid, 95% acetonitrile/water (v/v). Mobile phase B: 10 mM ammonium acetate, 0.1% acetic acid, 50% acetonitrile/water (v/v). Gradient: t = 0.0 min, 1% B, t = 1.0 min, 1.0% B, t = 3.0 min, 15.0% B, t = 6.0 min, 50.0% B, t = 9.0 min, 95% B, t = 10.0 min, 95% B, t = 10.5 min, 1% B, t = 15.0 min, 1% B. For both ionization modes, the column temperature was 35 °C and injection volume was 5 μL. Autosampler compartment was kept at 7 °C. Tune file settings: Sheath gas flow rate 30 (arbitrary units), aux gas flow rate 13 (arbitrary units), spray voltage 3.2 kV (2.7 kV negative ionization mode), capillary temperature 350 °C, aux gas heater temperature 400 °C.

Mass spectrometric settings

Full scan mass resolution 35,000, Scan range 70 to 1050 m/z, ACG target 1 e6, spectrum data type: profile. dd-MS2 discovery mode, MS2 mass resolution 17,500, isolation window 3.0 Da, collision energy 30 eV.

MS-data pre-processing and data analysis

Post-acquisition data processing was accomplished using Compound Discoverer™ (Thermo Scientific™) software 3.2, commonly used for untargeted metabolomics analyses. In brief, raw data was uploaded and processed; missing values imputation was automatically performed by the software. Then, features with more than 20% missing values across sampled were filtered.

The obtained mass to charge ratio (m/z), including MS/MS fragments and retention time were searched against large spectral libraries for metabolite annotation. After data evaluation and cleaning process, including background spectral filtering (noise elimination) and quality control (QC)-based filtering (metabolite features with CV ≥ 30% were removed), the dataset was suitable for statistical analysis. First, Principal Component Analysis (PCA), as an unsupervised statistical method, was utilized. Then, the screening for biomarkers was done by step-wise filtering of Adj P-value < 0.05, and different fold change ratios as detailed in the results section. Data exploration also included further investigation of the affected metabolic pathways by the online tool of MetaboAnalyst 5.0 software21. Finally, the detected potential biomarkers were confirmed manually by spectra evaluation and comparison to the fragmentation pattern and retention time of commercial standards or online library mzCloud™ and commercial standards by Sigma. Adjusted P-value < 0.05 was reported as significant.