Urine and Fecal ¹H-NMR Metabolomes Differ Significantly between Pre-Term and Full-Term Born Physically Fit Healthy Adult Males

2.5.1. Taxonomy and Microbial Diversity of Intestinal Tract

In contrast to the observed differences at the urine and fecal metabolomics levels described above, the taxonomic level of information did not result in significant differences between groups (PERMANOVA; p > 0.05; n permutations > 5000). In addition, based on the taxonomic data 181,020 JADBio models were trained using an extensive tuning effort, but no reliable biomarker or trained model could be obtained. In general, three different kingdoms were detected in all samples (archaea, bacteria, and DNA viruses). The average relative abundance of bacteria was lower in the preterm than in the full-term group (92.3% vs. 77.5%). The abundance of DNA viruses was higher in the preterm group (7.2% vs. 20%) (Figure S9a). The archaeal kingdom was least prevalent in both groups (0.5% in the control group vs. 2.5% in the preterm group). This lack of significant differences at taxonomic levels was previously attributed to large interpersonal differences between participants [58,59]; however, in this study variance within the full-term control group was at least two times larger than that observed in the preterm group and hence significantly higher (p < 0.05) based on the analysis of 3D coordinates after nmMDS and PCoA analysis. This points to the existence of an overarching effect shared by all preterm participants absent from the matching control full-term group.

We further point out that both groups contained matching groups of young healthy physically fit participants, in contrast to past studies exploring the differences between the preterm and general population [60]. Our results point to the fact that some microbiome-related characteristics within the preterm group were apparently shared to a larger extent within the preterm group in comparison to the full-term control group, signifying the existence of differences in the microbial makeup due to differences in the physiology of the host. The existence of an overarching effect shared by all preterm participants absent in the control group was evident from the significantly higher Shannon diversity (p < 0.05) (Figure S9b) in the preterm group including other diversity estimates that differed significantly between the two groups (Table S7). These two observations suggest the existence of an interplay between the increased similarity of the major preterm taxonomic categories and the diversity of a smaller highly variable list of taxa not well shared between the preterm participants in this study. It is easily envisioned that additional environmental factors shape the gut microbiome within the preterm group, making the preterm group a narrower subset of the otherwise healthy human gut.

2.5.2. Functional Genes of Human Gut Microbiome

Based on the 198,306 categories describing gene family data, in the preterm and control groups, 90,510 models were trained using an extensive tuning effort in search of biologically meaningful discriminative variables between the preterm and full-term control groups. The entire list of features was used to build and validate a trained model that achieved insignificant validation performance with an AUC and other metrics. Consequently, no significant differences could be identified between the two datasets at the level of functional gene lists.

2.5.3. Enzymatic Reactions Taking Place in Human Gut

For the aggregation of functional gene information into enzymatic reactions (Figure 4a), again 90,510 models were explored using an extensive tuning effort. The best model was ridge logistic regression with the penalty hyperparameter lambda = 1.0, with an area under the curve (AUC) value of 0.992. In addition to AUC, all other thresholds were also statistically significantly different from the baseline. Features were selected based on the Test-Budgeted Statistically Equivalent Signature (SES) algorithm with the following hyperparameters: maxK = 3, alpha = 0.1, and budget = 3 *n_vars. RXN-15378, RXN-14971, RXN-21393, and RXN-21394 were equally selected as the most important for discriminating between the preterm and control groups and were increased in the preterm group (Figure S10) (p < 0.05). All of the above reactions represent the enzymatic reaction succinate dehydrogenase based on the BioCyc website ([61], accessed on 23 January 2022). We used RXN-15378 to validate the trained model and obtained a validation performance with an AUC of 0.931. Succinate is a metabolite produced by both host and microbial cells and accumulates under conditions of inflammation and microbial imbalances in the intestinal tract [62]. Succinate was shown to accumulate in areas of inflammation and metabolic stress [63] and can have tissue specific but also systemic effects as a proinflammatory signaling molecule [62,64,65,66]. Although gut microbes represent the predominant source of succinate, it is typically rapidly consumed in the production of propionate, one of the major short chain fatty acids, by Bacteroides spp., Prevotella spp. and some members of Firmicutes [62,67]. Although the mucosal uptake of succinate as a charged molecule over the mucosal epithelia is significantly higher in the small intestine, it takes place to various extents throughout the length of the intestinal tract and requires sodium dependent transport proteins [62]. In addition to the internalized succinate provided by the microbiome, succinate also accumulates within cells under conditions of low oxygen as a metabolic signature of hypoxia, generating HIF-1α to regulate cellular responses and adapt to a low oxygen environment. At normoxia, HIF-1α is regulated by posttranslational hydroxylation and degradation by prolyl-hydroxylase activity that converts alpha-keto-glutarate to CO₂ and succinate while inactivating HIF-1α. Excess uptake of microbiome produced succinate results in higher levels of intracellular succinate that can slow down prolyl-hydroxylase activity through product inhibition and result in an additional activation and stabilization of HIF-1α beyond its response to hypoxia itself, which can significantly augment the LPS-induced expression of proinflammatory cytokines [62,68].

2.5.4. Metabolic Pathways Observed in Human Gut

Based on the metabolic pathway data (Figure 4b), 60,310 models were trained, with extensive tuning effort. The best model was ridge logistic regression with a penalty hyperparameter of 100 and an area under the curve (AUC) of 0.981. In addition to AUC, all other thresholds were also statistically significantly different from baseline. Features were selected on the basis of Lasso feature selection with a penalty = 1.5. On the basis of the MetaCyc website ([69,70,71]; accessed on 23 January 2022), the most important metabolic pathways were PWY-7456 (β-(1,4)-mannan degradation), PWY-7323 (superpathway of GDP-mannose-derived O-antigen building blocks biosynthesis), GLYCOLYSIS-TCA-GLYOX-BYPASS (a superpathway of glycolysis, pyruvate dehydrogenase, TCA, and glyoxylate bypass), P221-PWY (octane oxidation), and PWY-5173 (unclassified). These pathways were the most important for distinguishing the preterm group from the control group and were all increased in the preterm group (Figure S11). The entire set of selected features was used to validate the trained model and achieved a validation performance with an AUC of 1.00. The relative frequency of this response was significantly increased in the preterm group, which was also confirmed by the t-statistic (p < 0.05). A set of selected features was used to validate the trained model and achieved validation performance with an AUC of 1.00.

The β-(1,4)-mannan degradation (PWY-7456) belongs to Bacteroides fragilis in the human intestinal tract and is essential for mucosal integrity and host nutrition [72,73]. Degradation of mannan by either Bacteroides dorei or Fecalibacterium prausnitzii and Roseburia intestinalis promotes the growth of Lactobacillus helveticus and Bifidobacterium adolescentis, which have probiotic properties and promote the synthesis of short chain fatty acids [74] or promote the growth of commensal microbes [75,76].

The superpathway of GDP-mannose-derived O-antigen building blocks biosynthesis (PWY-7323) is involved in lipopolysaccharide (LPS) biosynthesis. Only gram-negative bacteria have LPS, and O-antigen is the part that extends the polysaccharide away from the cell surface and triggers the host cell immune response [77,78]. Gram-negative bacteria observed in preterm infants cause serious infections, such as sepsis [79] coupled with the absence of MD -2 (a protein responsible for the recognition of LPS), which leads to a higher risk of developing intestinal diseases in adults born preterm due to the impaired recognition of LPS in the past [78]. Elevated LPS levels may also contribute to inflammaging (chronic, low-grade inflammation that develops with age) [80]. This also fits our observation that microbially produced succinate coupled with hypoxia can significantly augment LPS-induced expression of the proinflammatory cytokines [68].

GLYCOLYSIS-TCA-GLYOX-BYPASS (a superpathway of glycolysis, pyruvate dehydrogenase, TCA, and glyoxylate bypass), is a superpathway that was significantly overrepresented in the preterm group. It integrates some of the fundamental components of energy metabolism, starting with a hexose sugar and ending with CO₂ and several forms of highly reducing metabolites that can be used for adenosine triphosphate (ATP) generation. Even though acetyl-CoA is shown in this superpathway as a product of the glycolysis pathway, it is also generated by the degradation of fats and proteins and by the fermentation of many metabolites. This superpathway includes the glyoxylate cycle, which bypasses those steps in the TCA cycle that lead to a loss of CO₂, and operates in bacteria. The increased energy production in the preterm microbiome apparently coincided with the general characteristics of the preterm individuals, such as increased oxidative stress, elevated antioxidant activity, and NO metabolism in acute exercise as described above [7]. It is possible to suggest that the intestinal conditions experienced by the gut microbiome exerted additional stress on the microbial functioning as well. Further research is needed to corroborate this notion.

In line with our observation of differences in microbiome functioning, P221-PWY (octane oxidation) was shown to increase with the Westernization of the human gut and lifestyle [81,82] as volatile organic compounds, including octane, were found either in the exhaled air or feces of human subjects with diverse medical conditions associated with oxidative stress and chronic inflammation, including lung cancer [82,83] obstructive sleep apnea [84], gastrointestinal diseases [85] and NAFLD [86]. Previous studies demonstrated that several alkane-degrading bacteria were capable of using diverse compounds as a carbon source in addition to alkanes [87], which are further oxidized to fatty acids via the bacterial β-oxidation pathway (BioCyc ID: P221-PWY). The key process in octane oxidation is the alkane hydroxylase system that introduces molecular oxygen in the C1 atom of the hydrocarbons at the expense of NADH to yield primary alcohols [88] that were further linked to liver associated diseases.

The acetyl-CoA biosynthesis (PWY-5173) pathway involved in carbohydrate metabolism was also significantly increased in the preterm group. The resulting acetyl-CoA acts as a precursor in the synthesis of intestinal short chain fatty acids including butyrate and acetate [89], that are important in maintaining gut health [90]. The increased levels of acetyl-CoA biosynthesis fit nicely with the other pathways observed in this study that either contribute or consume mass flow related to this reaction. As both preterm and full-term groups were composed of healthy young physically fit males differing significantly in acetyl-CoA biosynthesis, our findings support a recent report on this pathway being one of the most variable pathways in a survey of subgroups of elite Irish athletes [91].

These overall results of the metabolic pathway analysis point to the fact that (irrespective of the heterogeneous makeup of the underlying microbiome taxonomy within the individual participant) the complex coordinated adjustments to the metabolism of the microbiome nevertheless take place and can be robustly reproduced from the integration of the sequencing information as described in this study [91,92,93,94] and can be linked to physiologically meaningful differences between groups reported before [24].

2.5.5. Predicted Water- and Lipid-Soluble Intestinal Metabolites

Our last layer of information dealt with the extended analysis of sequencing data towards water- and lipid- soluble predicted metabolites utilizing relaxation-network analysis, which has been extensively trained and validated before [52,55]. In summary, 17 metabolites (out of 81) (Table S8) predicted with MelonnPan were detected also by ¹H-NMR in fecal samples, showing possible interaction between two systemic metabolisms (human and microbial). None of these metabolites were chosen by machine learning. Metabolites associated with the human gut microbiota (Figure 4c) were explored utilizing JADBio and 181,020 models were trained using extensive tuning efforts. The best model was a Support Vector Machine type C-SVC with a radial basis function kernel and hyper-parameter (cost = 10, gamma = 1.0), with an area under the curve (AUC) value of 0.976. In addition to AUC, all other thresholds were also statistically significantly different from the baseline. Feature selection was based on LASSO feature selection with a penalty = 0.25. Alpha-muricholate, putrescine, dimethyllysine, diacetylspermine, and C16 carnitine were significantly increased in the preterm group. In contrast, hydrocinnamic acid, fructose, glucose and galactose, chenodeoxycholate and deoxycholate were lower in the preterm group (Figure S12). When the trained model was applied to the test portion (30% of our total data set) validation performance with an AUC of 0.957 was obtained. In the following sections let us first review the predicted metabolites significantly increased in the preterm group.

Carnitine was increased in the preterm group and is associated with trimethylamine N-oxide (TMAO) production (Figure S12). TMAO is synthesized by the microbiota from trimethylamine (TMA), which in turn is formed from carnitine or choline. Increased choline content in the preterm group was also observed in urine metabolomics. These two molecules together (carnitine and choline), in conjunction with the microbiota, may be the most important cause of the increased likelihood of cardiovascular disease in the preterm group [8,95,96,97,98,99,100].

Putrescine and diacetylspermine are polyamines and important metabolites for the gut microbiota (Figure S12). Putrescine is synthesized by interspecies cooperation between Escherichia coli and Enterococcus faeacalis and is formed from arginine [101,102]. Elevated putrescine levels have been associated with an older gut microbiota [103], increased gut permeability, and elevated levels of inflammatory cytokines in mouse colon tissue [104]. Elevated putrescine levels led to activation of genes that regulate oxidative stress, which may lead to a parallel increased risk of developing metabolic syndrome [35,36,37,38] and irritable bowel syndrome [105].

Diacetylspermine as a polyamine metabolite was linked to cancer growth and its association with microbial biofilm formation. It is synthesized by bacterial acetylation and has been significantly upregulated in tissues with biofilms in animal models [106]. This suggests that microbial organization and biofilm formation capacity at the interface between the mucus layer and lumen might differ significantly between the preterm and full-term participants, an observation worth further exploration.

Dimethyllysine can be the end product of either host or microbial metabolism but currently little is known about its physiological roles for the host and microbiome in the Human Metabolome Database [107], ChemSpider (https://www.chemspider.com/ (accessed on 15 April 2022)) or FooDB (www.foodb.ca (accessed on 15 April 2022)) (Figure S12). A recent review of macronutrient metabolism by the human gut microbiome focusing on major fermentation byproducts and their impact on host health [108] reported that the major products of lysin were acetate, butyrate and cadaverine, hence linking this compound to the short- and long-chain fatty acid cycles associated with ulcerative colitis [109].

Alpha-muricholic acid was identified by the MelonnPan [55] relaxation network since its first use in the analyses of human samples analyzed using MelonnPan [110,111], suggesting a misclassification of rodent muricholic acid for cholic acid in humans in this approach. Nevertheless, irrespective of its MelonnPan supported assignment, it is evident that this secondary bile acid was identified at elevated levels in the preterm group, fitting into the framework of the distinct chemical makeup of the preterm gut in relation to fat metabolism and the metabolites reported in this study.

In addition to elevated metabolites identified by MelonnPan in the preterm group, the following metabolites were identified in significantly lower concentrations in the preterm group.

Deoxycholate (decreased in the preterm group) is another metabolite that interacts with microbes (Figure S12). Deoxycholate is a secondary bile acid. The human intestinal microbiota (Bacteroides intestinalis, Bacteroides fragilis, Escherichia coli) are involved in the production of secondary bile acids from primary bile acids, such as choline. Deoxycholate is also known to promote colon cancer. Because of the increased cholate levels in the preterm group, we would expect a greater likelihood of microbial metabolites associated with primary bile acids, as well as increased levels of the expected metabolites in the preterm group. In contrast, deoxycholate levels were decreased in the preterm group. This could also be due to the greater urinary excretion of cholate (cholate was increased in the urine of preterm infants). Bile acids also generally induce mitochondria ROS production. Preterm infants are challenged by ROS in the first few months of life, possibly implying that the systemic response is to increased urinary excretion of bile acids in preterm infants [112,113,114,115].

Hydrocinnamic acids are a major class of phenolic acids from dietary fiber with the characteristic phenylpropanoid C6-C3 backbone that were significantly decreased in the preterm group (Figure S12). Although the polyphenol–gut microbiota interactions and their impact on human health have been known for decades, there is great inter-individual variation caused by the different individual capabilities of processing, absorbing and using these compounds effectively [116]. In light of the physiological differences between the two groups analyzed in this study, it seems plausible that differences exist also in the extent of the utilization of these polyphenols in the preterm group. In addition, lower levels of hydrocinnamic acid were observed in patients with Crohn’s disease and ulcerative colitis compared to the healthy cohort. Lower levels of hydrocinnamic acid in the preterm group may lead to increased levels of circulating BCAAs, which in turn predisposes preterm born individuals to metabolic syndrome and cardiovascular disease [117,118,119,120].

The lower levels of reducing sugars fructose, glucose and galactose, in the preterm group corresponded with a greater capacity to produce short chain fatty acids (Figure S12). The metabolic reactions and predicted metabolites jointly suggest the existence of a larger metabolic flow-through of the preterm microbiome in comparison to the full-term group, pointing to significant differences in the environmental setup in the preterm gut.

In conclusion, the results presented here constitute the first report on the differences in the urine and fecal metabolomes between preterm and full-term groups of physically fit healthy young males. Clear differences were identified in the urine and fecal metabolomes next to the metabolic pathways, suggesting that systemic differences between the two groups affect the metabolism of the host as well as intestinal tract parameters and that of the underlying microbiome (Figure S13). One has to realize that studies with female participants are lacking and not many studies with sufficient statistical power were reported so far to close the gap. With the concomitant methodological developments presented in this study, the exploration of the more complex female metabolome and responses to inactivity and hypoxia can be commenced in a comparable way [35].

Spectra Processing

NMR spectra were preprocessed with an internal script and prepared for identification with the Chenomx Compound Library, extended to the Human Metabolome Database [41,107], giving access to the chemical shift profiles of 674 compounds used in the analyses. Chemical shifts of 647 compounds were used for the identification of metabolites observed in our study. The resulting spectra were then analyzed with targeted quantitative metabolomics using Chenomx NMR Suite version 8.6 (Chenomx, Inc., Edmonton, AB, Canada). ChenomX profiler was used for randomized spectral fitting. All spectra were processed in the same way by spectral deconvolution and once metabolites were identified, urine and fecal data matrices were established, assigning 0 to a particular metabolite not detected in all samples.

Sequence Processing

Paired reads obtained from Macrogen were analyzed using our in-house pipeline for metagenomics sequence processing—Metabakery (in preparation). Metabakery is a re-implementation of the BioBakery [52] workflow using (https://huttenhower.sph.harvard.edu/kneaddata/, accessed on 7 April 2022) for quality control, MetaPhlAn [53] for taxonomy analysis (bacteria, archaea, fungi, protozoa, and viruses), and HUMAn3 [53] for functional genes, enzymatic reactions, and metabolic pathways determination. Additionally, the MelonnPan was used for the prediction of metabolites. Metabakery is implanted as a singularity image and prepared on high computing performance clusters. The analyses running MetaBakery were performed on a dual Xeon system with 32 CPU cores (64 hyperthreads), 512 GB RAM and 6 TB SATA hard disc at the Faculty of Electrical Engineering, University of Ljubljana.

3.6.1. Statistics

First, the software PAST [129] was used for PERMANOVA. All obtained data matrices (NMR metabolomes—identified fecal and urinary metabolites at micromolar concentrations, microbial taxonomy, gene families, enzymatic reactions, metabolic pathways and predicted microbial metabolites) were analyzed in the same way. Each determined parameter was analyzed in three different ways as previously described [35,36,37,38,56]: (i) by dividing the measured concentration by the concentration of all metabolites in that sample; (ii) Box-Cox; or (iii) log(x + 1) transformed. The significance of the metabolic differences and microbial entities between the different sample groups were tested using ANOSIM, and NP-MANOVA, and expressed as the overlap in the non-metric multidimensional scaling (nm-MDS) trait space (using Euclidean distance measures). The stress function was used to select the dimensionality reduction, while Shepard’s plots were used to describe the correspondence between the target values and the obtained ranks. In addition, PCoA and PCA were performed on metagenomic data. Benjamini–Hochberg significance correction for multiple comparisons was used as previously described [130].

Second, for MetaboAnalyst [40], a log or cube root transformation was used in conjunction with mean or Pareto scaling as implemented in MetaboAnalyst, followed by supervised classification using the partial least squares discriminant analysis (PLSDA) method, random forest (RF), and pathway enrichment analysis. The PLSDA results were cross-validated with a caret package implemented in MeatboAnalyst. The major metabolites identified by PLSDA were determined according to the variable importance in projection (VIP). The randomForest package implemented in MetaboAnalyst was used for supervised classification between different groups of interest. The main features defined by RF were ordered according to the mean decrease in classification accuracy. Hierarchical clustering was performed according to the VIP scores to obtain a heat map representing the differences in metabolic profiles between samples and groups. Euclidean distance, Pearson’s correlation and Spearman’s correlation were used as similarity measures and Ward’s linkage was used as a clustering algorithm. MetaboAnalyst and gplot were used to generate graphs.

KEGG libraries for human metabolic pathways were used for metabolic pathway and enrichment analysis. For topological analysis, the globaltest analysis method and relative Betweenness centrality were used. Significant pathways were determined using the raw p-value, Holm–Bonferroni p-adjusted value, and adjusted p-value using the false discovery rate. The effect of pathways was calculated using the pathway topology analysis.

Metabolite Set Enrichment (MSEA) was used to identify biologically significant patterns between quantitative metabolome data from different groups. The names of compounds in Human Metabolome Database (HMDB) were used for linkage to the KEGG database. Enrichment analysis was performed using the globaltest package implemented in MetaboAnalyst. The enrichment ratio was calculated by dividing observed hits and expected hits.

3.6.2. JADBIO Auto Machine Learning

Just Add Data Bio (JADBIO), a web-based machine learning platform for analyzing potential biomarkers [57], was used to search for biomarkers. The JADBIO platform was developed for predictive modeling and providing high-quality predictive models for diagnostics using state-of-the-art statistical and machine learning methods. Personal analytic biases and methodological statistical errors were eliminated from the analysis by autonomously exploring different settings in the modeling steps, resulting in more convincing discovered features to distinguish between different groups. JADBIO with extensive tuning effort and six CPUs was used to model different dataset choices in addition to the features observed in samples of all groups from different projects by splitting the total data into a training set and a test set in a 70:30 ratio. The training set was used to train the model and the test set was used to evaluate the model [39,56].

To assess the classification of the model, a receiver-operating characteristic curve (ROC curve) was constructed for all studied groups, plotting the true-positive rate (sensitivity) against the false-positive rate (1-specificity). Individual conditional expectation plots (ICE) showed the nature of the contribution of each feature characteristic to the model. All obtained models can be run locally using a Java executor.