Introduction

During the transition period, the higher nutrient demands for fetal growth and milk production cause important metabolic adjustments to support energy requirements1,2. In parallel, a reduction in dry matter intake is established (DMI), leading to a negative energy balance (NEB)3. The main body reserves mobilized to support energy requirements and glucose synthesis are muscular and adipose tissues4. Nonesterified fatty acids (NEFA) are derived from the mobilization of adipose tissue and follow four different energy production pathways in the liver: (a) complete oxidation via the tricarboxylic acid cycle (TCA); (b) incomplete oxidation to generate ketone bodies and energy; (c) export from the liver in the form of very low-density lipoprotein (VLDL) or (d) esterification to triacylglycerols and accumulation within hepatocytes5,6.

However, the liver cannot metabolize the entire amount of NEFA derived from reserve deployment, generating an excess of ketone bodies (β-hydroxybutyrate (BHB), acetoacetate and acetone). This makes dairy cows more susceptible to metabolic diseases such as ketosis7. Ketosis is characterized by an increase in ketone bodies in the blood, urine and milk, which often occurs during early lactation8. The β-hydroxybutyrate concentration is used for diagnosis in high-yielding dairy cows, with a cutoff of blood BHB above 1.0–1.4 mmol/L for subclinical ketosis or hyperketonemia without clinical signs9,10,11. The presence of this disease is associated with economic losses due to the risk of abomasum dislocation, reproduction disorders, infectious disease, reduction in milk production and a higher risk of culling animals9,11.

Metabolomics is an analytical approach that aims to simultaneously measure the entire metabolite profile of a biologic sample12. The metabolites represent the final product of cellular processes in response to environmental changes13, so their evaluation could offer a close and complete view of how dairy cows react to metabolic stress in the transition period14. Among the analytical tools used for metabolomics investigations, proton nuclear magnetic resonance spectroscopy (1H-NMR) requires minimal sample preparation and provides robust and highly reproducible information12, which compensates for its lower sensitivity in comparison to other analytical platforms. Serum is a common biological fluid used in this type of analysis because it is easy to sample and minimizes animal welfare problems. Moreover, it can provide metabolites from all organs15 to associate them with metabolic disease16.

Progressive changes in the metabolomic status of animals during hyperketonemia are measurable. The aim of this study was to evaluate the serum metabolome analyzed using 1H-NMR in dairy cows with different levels of BHB.

Results

Main characteristics

The serum BHB and NEFA values presented significant differences among each of the groups analyzed (Table 1). In detail, the G1-G0 difference in BHB was 22% of the G2-G0 difference. Focusing on NEFA, Groups G0, G1 and G2 showed mean values of 0.23, 0.34 and 0.62 mEq/L, respectively, so that the G1-G0 difference was 35.9% of the G2-G0 difference.

Table 1 Least square means and standard error of the mean (SEM) of main characteristics for each group.
Full size table

Serum metabolome profile and robust principal component analysis (rPCA)

In the serum metabolome, fifty-seven metabolites were identified (Table 2). Among the identified metabolites, thirteen were significantly different between Groups G0 and G2 according to univariate analysis: glutamate, proline, serine, aspartate, isovalerate, and choline showed a significant reduction, whereas 3-hydroxybutyrate, 3-hydroxyisobutyrate, acetate, succinate, 2,3-butanediol, methanol, and methylsuccinate showed a significant increase. The calculated rPCA model is shown in Panels A–C of Fig. 1. The first principal component (PC1) of its scoreplot, accounting for 66.8% of all of the samples’ variance explained, nicely summarizes the overall differences between samples of Group G0, with low PC1 scores, and samples from Group G2, with high PC1 scores. The loading plot showed that the molecules mostly representative of the G0 group were choline, glutamate, proline, and aspartate, while molecules mostly representative of the G2 samples were 3-hydroxyisobutyrate, 2,3-butanediol, methylsuccinate and methanol. When samples pertaining to Group G1 were projected in this 12 molecules’ space (Fig. 1D,E), they appeared, along PC1, between Groups G0 and G2, with their distance from G0 being 36.3% of the G0–G2 distance.

Table 2 Mean values and standard error of the mean (SEM) of metabolite concentrations expressed in µmol/L and their p value corrected by the Bonferroni method.
Full size table
Figure 1
figure 1

rPCA model built on the space constituted by the concentration of the significant metabolites listed in Table 2 for Groups G0 and G2. In the scoreplot (A), samples from Groups G0 and G2 are represented with black squares and green triangles, respectively. The wide, empty circles represent the median of each sample group. The position of the samples along PC1 is summarized in the boxplot (B). The loading plot (C) reports the significant correlation between the concentration of each substance and its importance over PC 1 (p < 0.05). The G1 group is represented as red circles in scoreplot (D) and boxplot (E).

Full size image

In addition to the 13 significantly different metabolites, 11 seemed to show trends toward significance: lysine, alanine, arginine, formate, pyruvate, and dimethylsulfone were reduced, whereas isoleucine, valine, ethanol, trimethylamine-N-oxide (TMAO), and acetone were increased. To define the overall underlying structure, a rPCA model was built through the two-step procedure mentioned above (Supplementary Fig. S1). Samples G0 appeared at low scores along PC1, which represented 56.1% of all of the samples’ variance explained by the model, while samples G2 appeared toward positive values. The molecules mostly representative of the former were dimethyl sulfone, myo-inositol and pyruvate, while those mostly representative of sample G2 were valine and isoleucine. Samples G1, projected over this model, were superimposed with those from Group G0. Interestingly, while samples G2 were not different from samples G0 when the 11 molecules were observed singularly, the difference became significant when their overall trend was observed by rPCA.

Over representation analysis (ORA)

An ORA was used to gain information about the metabolic alterations that were possibly responsible for changes in the metabolome profile (Table 3). This highlighted that lipid metabolism was significantly (p = 0.016) altered by hyperketonemia. In fact, among the 7 molecules that could be referred to as the metabolism of lipids, 4 appeared to be significantly different between Groups G0 and G2. Studying lipid metabolism in deeper detail, the synthesis of phosphatidylserine seemed to be peculiarly overrepresented, with each of the two molecules observed (choline and serine) significantly altered. Another metabolic pathway influenced by hyperketonemia was glycosaminoglycan metabolism.

Table 3 Overrepresentation analysis (ORA) based on the metabolites significantly different between healthy (G0) and hyperketonemic (G2) animals.
Full size table

Discussion

During the transition period, cows develop NEB status due to a rapid increase in energy demand after the partum, high consumption of glucose by the mammary gland and a corresponding reduction of its availability, possibly causing hypoglycemia17. In response to this condition, cows undergo an increase in gluconeogenesis18.

The main glucogenic substrates are propionic acid, lactic acid, glycerol and glucogenic amino acids14. These substrates are used because sugars and starches are metabolized during ruminal fermentation, while they are common sources of glucose in monogastric species19. A deficiency of glucogenic substrates is an important risk factor for the pathogenesis of ketosis14. This disease often occurs between the second and seventh weeks postpartum in dairy cows, and subclinical ketosis or hyperketonemia is frequently observed20. A better knowledge of healthy and pathologic animals’ metabolome profiles could be useful to understand metabolic alterations associated with BHB increments9,21.

In agreement with previous investigations10, samples from Group G2 showed NEFA values significantly higher than those from Group G0, confirming the excessive lipid mobilization in bovines during peripartum. An excess of circulating NEFA can induce BHB production through incomplete fatty acid oxidation5,7. Oxidation of NEFA may lead to an increase in reactive oxygen species (ROS) in mitochondria due to the respiratory chain. Therefore, the greater concentration of NEFA and BHB may suggest a state of oxidative stress. In addition, NEFA may positively influence ROS production in neutrophils further increasing oxidative stress and influencing the immune response22,23,24. More interestingly, Group G1 showed NEFA values between those of G0 and G2, which were significantly different from both groups. This confirms that cows identified as healthy but with an increment of ketone bodies are indeed switching to a diseased status; consequently, they would need proper management and to be managed separately from healthy animals. However, studies focusing on better management of animals enrolled in the G1 group have not yet been conducted to authors’ knowledge.

In agreement with the literature11, in our study, the mobilization of body reserves was represented by an increase in the ketone bodies BHB and acetone in the G2 group, even if only the former was significant, due to a high variance in the latter. BHB is a product derived from muscular amino acid degradation or it can be synthesized from acetyl-CoA. It can be used as an energy source for other tissues, such as the brain and heart11. According to Zhang et al. (2017), it can influence the inflammatory response by reducing macrophage function and the production of proinflammatory cytokines.

Ketone bodies can be synthesized from ketogenic amino acids (lysine and leucine) and glucogenic and ketogenic amino acids (isoleucine)25,26. However, while isoleucine tended to increase from G0 to G2, lysine tended to show the opposite trend, in agreement with Sun et al. (2014) and Wang et al. (2016). These findings may suggest that lysine could be one of the first amino acids utilized as an energy supply. In addition, in humans, high levels of branched-chain amino acids (BCAAs; isoleucine, leucine, and valine) have been found to be connected to obesity and diabetes, with reports that an increase in BCAAs can lead to insulin resistance and, on the other hand, that higher insulin concentrations can promote protein synthesis27. In our study, the increase in BCAA concentrations (isoleucine and valine) in hyperketonemic cows could suggest that hyperketonemic cows may develop progressive insulin resistance proportionally to the BHB concentration. In a recent study on dairy cows, a single-dose duodenal infusion of leucine was not effective in stimulating an insulin response28. However, this study examined only six cows in late lactation and evaluated only the response to leucine, which was not significant in our study. Further investigation may therefore be necessary.

Methanol, which was significantly more concentrated in the G2 group, could be derived from methionine metabolism26 or from methane during the methane metabolism pathway. However, methionine does not exhibit concentration changes that would suggest an alteration in its metabolism. Methane is produced during ruminal microbial fermentation and is positively related to the mean retention time of particulates and liquid in the rumen29. The significant increase that we observed for methanol allows us to speculate that methane was more concentrated in the rumen, produced from other metabolites. As a first confirmation, Yanibada et al. (2020) found in the plasma of dairy cows that dimethylsulfone and formate are inversely related to methane ruminal production. A second confirmation comes from acetate, a fatty acid produced during ruminal fermentation30, which was significantly increased in Group G2. This fatty acid is positively related to methanogenesis31 and can be used by animals for energetic production when bound to coenzyme A32. Furthermore, we should consider that all groups were fed the same total mixed ration (TMR) with no differences in the diet administered.

The increase in methanol can also be related to the increase in ethanol concentration because ethanol can inhibit methanol utilization by microorganisms33. Ethanol is a product of anaerobic fermentation in the rumen by yeast and bacteria34 that can act as an agonist of the GABA receptor with a consequential inhibitory effect. This metabolite can be used by ruminants for energetic production through its conversion to acetate35.

The serum concentration of 2,3-butanediol, significantly higher in Group G2, agrees with the observations by36, who described this metabolite as produced from 2-butanone through acetoin in several ketoacidotic conditions. In humans, 3-hydroxyisobutyrate serum levels were also found to be directly correlated with the level of ketoacidosis37. Experiments on rats suggested that under ketoacidotic conditions, this molecule is produced from valine37. Overall, the increment of methanol, acetate, ethanol, 2,3-butanediol, and 3-hydroxyisobutyrate and the reduction of dimethylsulfone and formate seem to suggest a pathologic increase in ruminal fermentation in hyperketonemic cows.

Choline may be consumed along the metabolic pathway toward glycine synthesis and converted to betaine and then to dimethylglycine. In our study, we observed all of these metabolites, but only choline appeared to be significantly altered by the levels of BHB or hyperketonemia. A more likely fate of choline is therefore its conversion to TMAO, which in our study showed a trend toward an increase and a connection with hyperketonemia. According to Xu et al. (2016), TMAO is a marker of oxidative stress because it reduces glycolysis and enhances β-oxidation of fatty acids. As an alternative involvement of choline in lipid metabolism, Sun et al. (2014) noticed that choline supports fatty acid transport and reduces the risk of hepatic lipidosis. The significant reduction in choline may therefore suggest an alteration of lipid transport with an enhanced risk of developing hepatic lipidosis.

The above pieces of information taken together (NEFA, BHB, choline, and TMAO) could reinforce the idea of oxidative stress due to an enhancement of β-oxidation of fatty acids, and it could, in turn, explain the high risk of hepatic lipidosis in the hyperketonemic group.

Methylsuccinate, significantly increased in Group G2, has been previously related to ketotic conditions38, probably in relation to acyl-CoA dehydrogenase activity in the β-oxidation of fatty acids39. In accordance with the studies of40 and41, high values of methylsuccinate in urine can be related to genetic disorders such as acyl-CoA dehydrogenase deficiencies. Methylsuccinate is a metabolite related to succinate, an intermediate of TCA, also significantly increased in Group G2.

In addition, isovaleryl-CoA dehydrogenase is a branched-chain dehydrogenase42, which degrades isovaleryl-CoA during the leucine cycle. In the next step, the biotin-dependent enzyme methylcrotonyl-CoA carboxylase (MCC) completes the degradation of isovaleryl-CoA to acetoacetate and acetyl-CoA43. Furthermore, other biotin-independent pathways may be involved in leucine catabolism, particularly MCC, which is not available, but they require high energy43. Isovaleryl-CoA can be converted into isovalerate, a branched-chain saturated fatty acid. In our study, the isovalerate concentration was significantly decreased in the G2 group, with a trend opposite to that of acetate, as reported recently44.

This could suggest that isovaleryl-CoA was metabolized by other pathways that can be biotin-independent, with a loss of energy. As previously reported by45, the use of biotin can improve dry matter intake, milk yield, and body energy metabolism due to the gluconeogenesis pathway. Furthermore, biotin supplementation may be useful to improve isovaleryl-CoA metabolism for energy production.

Alanine and serine are glucogenic amino acid substrates for pyruvate synthesis46. In our study, alanine showed a tendency toward significance, and the change in the serine concentration was significant. The tendency of these metabolites was a progressive decrease with increasing BHB concentration. Furthermore, serine is an important regulator of glutathione synthesis and it is involved in oxidative stress management47. The reduction in serine may suggest an inability of the cows to handle a possible oxidative stress state in the early lactation period.

Pyruvate and its precursors alanine and serine showed an inverse relationship with the BHB concentration. Alanine represents one of the major resources for gluconeogenesis48. Overall, this seems to suggest a progressive involvement of pyruvate in gluconeogenesis, with a reduction in its concentration, due to the lack of different substrates in Group G2 compared to Group G0.

The first intermediate of TCA is oxaloacetate, which can be used for pyruvate synthesis. Aspartate is a precursor of oxaloacetate46. Although oxaloacetate was not identified in our study, aspartate showed a significant reduction in the G2 group. This finding may suggest a possible reduction in oxaloacetate in the hyperketonemic group due to the decrease in precursors involved in its synthesis.

Glutamate, proline, and arginine are glucogenic amino acids leading to α-ketoglutarate synthesis46. All of these amino acids are related to glutamate production49. Furthermore, proline and arginine are linked in the metabolic pathway of arginine and proline metabolism. In our study, proline and arginine concentrations decreased in G2. The progressive reduction in proline concentration was in agreement with the study of Wang et al.16 on ketotic cows. Proline and its products are components of collagen biosynthesis, which contributes to structure, strength and tissue integrity50. Furthermore, proline has been found to act as a weak agonist of the glycine receptor and both N-methyl-d-aspartate (NMDA) and non-NMDA ionotropic glutamate receptors. Arginine is a urea cycle amino acid51 that participates in nerve signal transduction, and its reduction is related to ketosis17. The reduced concentration of these metabolites during the glucogenic process could suggest an influence of the urea cycle, signal transduction and recovery from the inflammatory process in hyperketonemic cows.

Isoleucine and valine can be used to synthesize succinyl-CoA and succinate, two other intermediates of TCA46. As previously mentioned, succinate showed a significant increase in G2. Higher concentrations of succinyl-CoA led to a block of citrate synthase and α-ketoglutarate dehydrogenase52. The following step is the oxidation of succinate to fumarate through succinate dehydrogenase53. This context could lead to a reduction in succinyl-CoA production, an increase in succinate synthesis and its utilization for TCA. However, fumarate concentrations did not show significant changes in the G2 group, as did the succinate concentration, which was fumarate’s precursor in TCA. This result suggests that there was a block in the oxidation of succinate. Succinate dehydrogenase is an enzyme involved in TCA and in the electron transport chain, which is comprised of iron-sulfur protein and flavoprotein (FAD) subunits52,53. A block of this enzyme may be due to reduced levels of minerals or altered concentrations of FAD. Furthermore, acyl-CoA-dehydrogenases involved in fatty acid oxidation contain FAD39. Further studies specifically focused on this step of TCA may be advantageous to analyze this aspect.

The overrepresentation analysis, applied to the metabolites overall identified and to those significantly different between Groups G0 and G2, highlighted lipid metabolism as significantly (p = 0.016) altered by ketosis. In fact, among the 11 molecules that could be involved in the metabolism of lipids, 4 appeared to be significantly different between Groups G0 and G2. Studying lipid metabolism in greater detail, the synthesis of phosphatidylserine seemed to be particularly represented, with each of the two molecules observed to be significantly altered. Phosphatidylserine is synthesized by facilitating the exchange of l-serine with the choline head group in phosphatidylcholine and with the ethanolamine head group in phosphatidylethanolamine.

Furthermore, glycosaminoglycan metabolism was influenced in the hyperketonemic group. Glycosaminoglycans (GAGs) are long, unbranched polysaccharides containing a repeated disaccharide unit comprised of a hexosamine. GAGs are located primarily in the extracellular matrix (ECM) and on cell membranes. GAGs participate in many important signaling events, such as neuronal growth, inflammation and development54.

Conclusion

Metabolomic analysis through 1H-NMR is a useful tool to achieve knowledge about metabolic profiling related to serum β-hydroxybutyrate modifications during the transition period in dairy cows. The metabolic state of our hyperketonemic cows suggests (1) a mobilization of body resources; (2) increased anaerobic fermentation; (3) alteration of lipid metabolism; and (4) a potential oxidative stress state. Furthermore, this metabolic profiling proposes the lack of glucogenic substrates and potential alteration of the electron transport chain are involved in ketosis in dairy cows. These findings indicate a possible alteration of inflammatory and healing processes. This study demonstrates that the metabolomic approach can be considered a significant means to achieve knowledge about dairy cow diseases and their pathogenesis.

Methods

The Ethics Statement was approved by the Animal Care and Use Committee of the University of Padua (ID number 91/2019—“BovineOmics” Projects). Animal care and procedures were conducted in accordance with the Guide for the Care and Use of Laboratory Animals and Directive 2010/63/EU for animal experiments (National law: D.L. 26/2014). This study was carried out in compliance with the ARRIVE guidelines. Informed consent was obtained from the owners for handling the animals and for the clinical activity of the Veterinary Teaching Hospital, University of Padua.

Animals

Forty-nine Holstein Friesian dairy cows between 15 and 30 days in milk were enrolled from a single high-yielding dairy farm located in the province of Padua (Italy). These animals were selected among dairy cows used for the experimental design of Fiore et al.10. The farm had 600 dairy cows and 400 lactating animals, with a mean production of approximately 10.000 kg/cow/lactation. All animals had a dry period of approximately 60 days, with no streaming-up. The same TMR was used for all enrolled animals (Table 4).

Table 4 Total mixed ratio (TMR) and feedstuffs used for all groups (G0, G1, and G2).
Full size table

A clinical visit was performed for each animal by veterinarians at the University of Padua. Animals with signs of clinical disease (metritis, retained placenta, abomasal displacement, milk fever, lameness, or mastitis) on the day of sampling or between parturition and the day of sampling were not accepted for this study. Data about parity, days in milk (DIM), milk yield and body condition score (BCS) with a scale of 1 to 5 points55 were recorded. The animals did not show a BCS > 3.5 at parturition or on the day of sampling.

Experimental design

A cross-sectional experimental design was used. Blood sampling was carried out in the late morning on the day of enrollment during the clinical examination. Samples were collected from the coccygeal vein with a vacutainer system for each enrolled cow. The samples were stored in Venosafe tubes (9 mL; Terumo Venosafe, Leuvel, Belgium) containing Clot Activator.

The Venosafe tubes were refrigerated at 4 °C and transported in a portable freezer (CoolFreeze CFX65 W professional, Dometic, Stockholm, Sweden—minimum temperature −22 °C) at the same constant temperature to the laboratory of the Department of Animal Medicine, Production and Health (MAPS) of the University of Padua (Italy) within 1 h of blood sampling. The samples were centrifuged at 3000 rpm for 10 min (Heraeus Labofouge 400, Thermo Scientific, Milan, Italy). Two aliquots of serum were obtained from each cow and were placed in 1.5 mL Eppendorf tubes. One aliquot was stored at −20 °C for biochemical analysis and the other was stored at −80 °C for metabolomic evaluations.

Blood analysis and group division

The biochemical analysis was performed using an automatic clinical chemistry analyzer (BT3500 Biotecnica instruments SPA, Rome, Italy). The serum BHB concentration was measured using β-hydroxybutyrate enzymatic kinetics (Randox, Milan, Italy; BHB, mmol/L); NEFA were assessed through the NEFA RX Monza test colorimetric method (Randox, Milan, Italy; NEFA, mEq/L). The glucose concentration was measured using a colorimetric method (Biotecnica Instruments SPA, Rome, Italy).

According to the serum BHB concentration, the animals were divided into three groups (Table 1): Group 0 (G0; 12 healthy animals with 4 primiparous, 2 secondiparous, and 6 pluriparous; BHB ≤ 0.50 mmol/L); Group 1 (G1; 19 healthy animals with 5 primiparous, 3 secondiparous, and 11 pluriparous; 0.51 ≤ BHB < 1.0 mmol/L); and Group 2 (G2; 18 hyperketonemic animals with 2 primiparous, 2 secondiparous, and 14 pluriparous; BHB ≥ 1.0 mmol/L). The BHB in the G2 group ranged from 1.00 mmol/L as the minimum value to 2.77 mmol/L as the maximum value.

Metabolomic analysis

The metabolomics investigation was carried out through an NMR analysis solution with 10 mM 3-(trimethylsilyl)-propionic-2,2,3,3-d4 acid sodium salt (TSP) in D2O set at pH 7.00 ± 0.02 by means of 1 M phosphate buffer containing 2 mM NaN3. TSP was used as an NMR chemical-shift reference, while NaN3 avoided microbial proliferation, as suggested by56.

Serum samples were prepared for 1H-NMR by thawing and centrifuging 1 mL of each sample for 15 min at 18,630g and 4 °C. The supernatant (700 μL) was added to 100 μL of NMR analysis solution. Finally, each sample was centrifuged as previously mentioned.

1H-NMR spectra were recorded at 298 K with an AVANCE III spectrometer (Bruker, Milan, Italy) operating at a frequency of 600.13 MHz, equipped with the software Topspin 3.5. According to56, the signals from broad resonances originating from large molecules were suppressed by a CPMG filter comprised of 400 echoes with a τ of 400 μs and a 180° pulse of 24 μs for a total filter of 330 ms. The water residual signal was suppressed by means of a presaturation technique. This setting employed the cpmgpr1d sequence, part of the standard pulse sequence library. Each spectrum was acquired by summing 256 transients using 32 K data points over a 7184 Hz spectral window, with an acquisition time of 2.28 s and a relaxation delay of 5 s.

The spectral phase was manually adjusted in Topspin, while the subsequent adjustments were performed in R computational language by means of a script developed in-house57. After the removal of the residual water signal, the1H-NMR spectra were baseline-corrected by means of peak detection, according to the “rolling ball” principle58, implemented in the baseline R package59. The signals were assigned by comparing their chemical shift and multiplicity with the Chenomx software library (Chenomx Inc., Canada, ver. 8.3).

The molecules of the first serum sample analyzed were quantified by means of an external standard by taking advantage of the principle of reciprocity60. Differences in water content among samples were then taken into consideration by probabilistic quotient normalization61. Molecule quantification was performed by means of rectangular integration, considering one of the corresponding signals free from interferences.

Statistical analysis

Statistical analysis was conducted with R software ver. 4.0.3 computational language57. Differences among the main characteristics (NEFA, glucose, BCS, DIM, parity, and milk yield) of the groups were tested through one-way ANOVA using the group as a fixed factor (p value < 0.05).

The concentration of each serum metabolite not normally distributed was normalized by Box and Cox transformation62. An initial comparison assessing the interactions of parity and DIMs versus groups was performed with a linear mixed model. However, no significant results were found, so it was decided to use only group as a fixed effect. Comparisons between the extreme groups (G0 vs. G2) were then made by t-tests, considering p=0.05 as the limit of significance. A post hoc pairwise comparison among metabolite concentrations was performed using Bonferroni correction. A t-test was also used to evaluate differences between the G0-G1 and G1-G2 groups. A trend for significance was considered for metabolites with 0.05 ≤ p value ≤ 0.10.

Robust principal component analysis (rPCA)63 was employed on centered and scaled data of significant metabolites as a means to summarize the structure of the data. rPCA was performed through the PcaHubert algorithm, implemented in the “rrcov” package. First, the algorithm detects outlying samples by computing their distance from the others along and orthogonally to the PCA plane. The optimal number of principal components (PCs) is finally determined. The rPCA model is summarized by a score plot and a correlation plot. The score plot shows the overall structure of the data by showing the samples in the PC space. The second plot shows the molecules that mostly determine the structure of the data by reporting the correlations between the concentration of each molecule and the PCs.

Through the MetaboAnalyst (ver. 5.0) web-based application64, the databases of PubChem (https://pubchem.ncbi.nlm.nih.gov/), Human Metabolome Database (HMDB; https://hmdb.ca/metabolites/) and Kyoto Encyclopedia of Genes and Genomes (KEGG; https://www.genome.jp/kegg/) were consulted to obtain a functional interpretation of the metabolites.

Metabolic pathways overrepresentation analysis (ORA) was performed by Fisher’s exact test and by employing Reactome (https://reactome.org) as a pathways’ database.