Glucuronidation Metabolomic Fingerprinting to Map Host-Microbe Metabolism

Abstract

Glucuronidation is an important detoxification pathway that operates in balance with gastrointestinal microbial β-glucuronidase (GUS) enzymes that can regenerate active metabolites from their glucuronidated forms. Although significant progress has been made in characterizing GUS enzymes, methods to comprehensively define the glucuronidome - the collection of glucuronidated metabolites - remain limited. In this study we employed pattern-filtering data science approaches alongside untargeted LC-MS/MS metabolomics to map the glucuronidome in urine, serum, and colon/fecal samples from gnotobiotic and conventional mice. Our findings reveal microbiome-driven shifts in the glucuronidome, highlighting how differential GUS activity can influence host metabolite profiles. Reverse metabolomics of known glucuronidated chemicals and glucuronidation pattern filtering searches in public metabolomics datasets exposed the diversity of glucuronidated metabolites in human and mouse ecosystems. In summary, we present a new glucuronidation fingerprint resource that provides broader access to and analysis of the glucuronidome. By systematically capturing glucuronidation patterns, this resource enhances unknown metabolite annotation efforts and provides new insights into the dynamic relationship between the host and bacterial biotransformation activities.

Figure 1

Detecting and identifying glucuronidated features from untargeted LC-MS/MS. (A) Analysis workflow. The dashed arrows connect the steps for creating in-house LC-MS/MS libraries of glucuronidated, phytochemical, and gut microbial metabolites. The solid arrows indicate the complementary workflow for detecting and identifying glucuronidated features from biological samples. The MS/MS libraries generated from the left side and public MS/MS libraries are integrated into the identification pathway starting at the alignment step. The dotted arrow indicates the integration point of the two workflows. (B) Example and MS/MS visualization of a MassQL query for identifying glucuronidated features in negative mode MS/MS data. Font colors in the query correspond to the MS/MS feature colors. (C) Molecular network of glucuronidated and phytochemical standards. Larger nodes were returned by MassQL as having MS/MS features of glucuronidation. Nodes are colored by common aglycone. Gray nodes do not share a common aglycone.

Figure 2

Acetaminophen Effects on Glucuronidome of Mice. (A) Principal coordinate analysis of all features with MS/MS spectra acquired in negative ionization mode in urine from acetaminophen (N-acetyl-para-aminophenol, APAP) and control treated mice (n=8) (B) Principle coordinate analysis of all urine glucuronidated features, defined as having a neutral loss with m/z=176.0321 or 194.0425 and fragments at m/z=113.0244 and 85.0295, in the urine of APAP and control treated mice (n=8). (C) Differential abundance of glucuronidated features from the urine of APAP compared to control treated mice. red=significant increase in abundance with APAP. blue=significant decrease in abundance with APAP. (D) Molecular network of APAP metabolites (negative ionization mode). Nodes are labeled with metabolite m/z. Edges connecting phase 2 metabolites are labeled with the mass difference between the two nodes. (E) MS/MS mirror plot of APAP glucuronide (top) and APAP (bottom) detected in mouse urine. Green features are common to both MS/MS spectra. The gray fragment at m/z 326.0882 is the precursor m/z for APAP glucuronide. MS/MS features used in the MassQL query for glucuronidated features are indicated with “@”. (F) MS/MS mirror plot comparing APAP glucuronide detected in mouse urine (top) to an analytic standard (bottom). Graph and network metabolite labels: a - APAP glcA; b - thiomethyl APAP glcA; c - APAP glcA; d - methoxy APAP glcA; e - thio APAP glcA; f - APAP sulfate; g - APAP; h - thio APAP; i - cystein-S-yl APAP; j - cystein-S-yl APAP [2M-H]⁻; k - APAP mercapturate glcA; m- APAP mercapturate; n - methyl-3-thioAPAP; o - methyl-3-thioAPAP sulphoxide; p - methyl-3-thioAPAP sulfoxide; s - S-methyl-3-thioAPAP. PCoA R² values were calculated by adonis tests. GlcA count differences were tested with poisson regression. Individual GlcAs were tested using Student’s T-test and considered significant with a Benjamini & Hochberg FDR<0.1 and log2FC>1.0 or log2FC<1.0. Volcano plots: blue=less abundant in FMT, red=more abundant in FMT, compared to germ-free, yellow=select annotations.

Figure 3

Conventional FMT shifts the glucuronidome from the colon to the serum and urine. (A) Principal Coordinates Analysis of bray curtis distances of colonic GlcA feature areas. (B) Number of colonic GlcAs detected with areas greater than noise threshold. (C) Volcano plot of the log2(fold change) of colonic GlcAs in FMT compared to GF. (D) Principal Coordinates Analysis of Bray Curtis distances of serum GlcA areas. (E) Number of serum GlcAs detected with areas greater than noise threshold. (F) Volcano plot of the log2(fold change) of serum GlcAs in FMT compared to GF. (G) Principal Coordinates Analysis of bray curtis distances of day 0 and day 7 urine GlcA areas. (H) Number of day 0 and day 7 urine GlcAs detected with areas greater than noise threshold. (I) Volcano plot of the log2(fold change) of day 7 urine GlcAs in FMT compared to GF. Metabolite names for panels C, F and I: a - daidzein 7-O-✉-D-glucuronide; b - ✉-muricholic acid glucuronide conjugate 4; c - genistein 7-O-✉-D-glucuronide; d - genistein 4’-O-✉-D-glucuronide; e - 3,5-dihydroxyphenylpropanoic acid 3-O-✉-D-glucuronide; f - naringenin-7-O-β-D-glucuronide; g - indoxyl β-D-glucuronide; h - phenyl β-D-glucuronide; i - p-cresol glucuronide; j - R,S equol 7-β-D-glucuronide; k* - dihydrogenistein glucuronide (* indicates tentative identification from unconjugated dihydrogenistein standard); m* - equol 4’-β-D-glucuronide (* indicates tentative identification, retention time shift from R,S equol 7-β-D-glucuronide). PCoA R² values were calculated by adonis tests. GlcA count differences were tested with poisson regression. Individual GlcAs were tested using Student’s T-test and considered significant with a Benjamini & Hochberg FDR<0.1 and log2FC>1.0 or log2FC<1.0. Volcano plots: blue=less abundant in FMT, red=more abundant in FMT, compared to germ-free, yellow=select annotations.

Figure 4

Low dose oral antibiotics shift the urine, colonic and serum glucuronidomes. (A-C) Counts of glucuronidated features detected in urine, colonic contents and serum, respectively, in control, gentamicin and vancomycin treated mice. (D) PCoA of urine glucuronidome after five days of vehicle, gentamicin or vancomycin treatment. (E-F) Volcano plot of the Log2 fold change (Log2(FC)) vs. -log(p-value) of urine glucuronidated feature areas after 5 days of gentamicin or vancomycin treatment, respectively, compared to control. (G) PCoA of colonic content glucuronidome after five days of vehicle, gentamicin or vancomycin treatment. (H-I) Volcano plot of the Log2 fold change (Log2(FC)) vs. -log(p-value) of colonic content glucuronidated feature areas after 5 days of gentamicin or vancomycin treatment, respectively, compared to control. (J) PCoA of serum glucuronidome after five days of vehicle, gentamicin or vancomycin treatment. (K-L) Volcano plot of the Log2 fold change (Log2(FC)) vs. -log(p-value) of serum glucuronidated feature areas after 5 days of gentamicin or vancomycin treatment, respectively, compared to control. Metabolite names for panels E, F, and L: a* - dihydrogenistein glucuronide (* indicates tentative identification to unconjugated dihydrogenistein standard) ; b - dihydro caffeic acid-3-O-β-D-glucuronide; c - R,S equol 7-β-D-glucuronide; d - dihydro ferulic acid 4-O-β-D-glucuronide; e - apigenin 7-glucuronide; f - phenyl β-D-glucuronide; g - kaempferol-3-glucuronide; h - indoxyl β-D-glucuronide. PCoA R² values were calculated by adonis tests. GlcA count differences were tested with poisson regression. Individual GlcAs were tested using Dunnett’s Test and considered significant with a Benjamini & Hochberg FDR<0.1 and log2FC>1.0 or log2FC<1.0. Volcano plots are treatment compared to control: blue=less abundant in treatment, red=more abundant in treatment, yellow=selected annotations.

Figure 5

Annotating Glucuronidated Features Through Molecular Networking. (A) Classical molecular networking subcluster of negative ESI heme-metabolite features. The structures are bilirubin glucuronide (left) and biliverdin glucuronide (right), with the glucuronic acid highlighted in yellow. (B) Classical molecular networking subcluster of positive ESI heme-metabolite features. (C) Classical molecular networking subcluster of soy-related flavonoid metabolites. The MassQL results for the glucuronidation pattern are shown as diamond shaped nodes. Putative annotations were achieved by MS/MS comparison to the connected non-glucuronide features for shared fragmentation. Two examples of aglycone:glucuronide MS/MS comparisons are included for biochanin A and formononetin. * indicates the putative annotation was later confirmed with standards. Panel C metabolite annotations: a - daidzein 7-β-D-glucuronide; b - daidzein; c - genistein; d - genistein sulfate; e - biochanin A sulfate; f - biochanin A; g - biochanin A glucuronide; h - daidzein 4’-β-D-glucuronide; i - genistein 7-β-D-glucuronide; j - hispidulin glucuronide (aka 6-O-methyl scutellarin); k - formononetin sulfate; m- formononetin; n - formononetin glucuronide

Figure 6

Microbial GUS driven glucuronidated and aglycone metabolite shifts. (A) Log2 transformed peak areas of colonic dihydrogenistein (aglycone) and dihydrogenistein glucuronide and urine dihydrogenistein glucuronide in germ-free (GF) mice and mice 7 days post fecal microbiota transfer (FMT). Hypothesis testing was performed by Student’s T-test. (B) Correlation plot of Schaedlerella_MGG38568_4_CTD GUS log transformed, normalized read counts and colonic aglycone dihydrogenistein log2 transformed peak areas, with linear regression fit line (y=2.7451x+22.1757). (C) Correlation plot of Schaedlerella_MGG38568_4_CTD GUS log transformed, normalized read counts and urine dihydrogenistein glucuronide log2 transformed peak areas, with linear regression fit line (y=2.2304x+25.5447). Adjusted R² and p-value generated using lm in R. (D -E) Volcano plot of log fold change (FC) by -log(p-value) for GUS genes in Gentamicin (gent) versus untreated control and vancomycin (vanc) versus control mice, respectively. The number in parentheses in the lower corners of each volcano plot indicate the number of positively and negatively differentially expressed genes. The point shapes indicate the loop class of each gene and color represents the statistical significance: blue=less abundant in treatment compared to control, gray=no statistical difference compared to control and red=more abundant in treatment compared to control. Hypothesis testing was performed with linear mixed effects models in R, using lmer for the main model and emmeans for pairwise comparisons between groups. Lower case letters indicate up-regulated CTD and Loop 1 GUS genes: a - Schaedlerella_MGG38568_3_CTD, b - Schaedlerella_MGG38568_1_CTD, c - Schaedlerella_MGG38568_4_CTD, d - X1XD42_69_sp011959925_Loop_1, e - Roseburia_MGG22730_3_CTD, f - Roseburia_MGG22730_4_CTD, g - Roseburia_MGG22730_5_Loop_1, h - Roseburia_MGG22730_1_CTD, i - UBA9475_MGG43629_CTD, j - Oscillospirales_unclassified_CTD. (F) Proportional peak area of naringenin (Nar) to naringenin-7-O-β-D-glucuronide (Nar-7-GlcA) after 30 min incubation of 0.2nM of 24 representative purified GUS enzymes and 50μM naringenin-7-O-β-D-glucuronide. GUS genes are clustered by loop classification. CTD - C-terminal domain, L1 - Loop 1, L2 - Loop 2, mL1 - mini-Loop 1, mL1,2 - mini-Loop 1,2, mL2 - mini-Loop 2, NL - No Loop and NTL – N-terminal domain. (G) Proportional peak area of naringenin (Nar) to naringenin-7-O-β-D-glucuronide (Nar-7-GlcA) after 24h incubation with Escherichia coli (E. coli) DSM 18039, which naturally expresses the L1 EcGUS used in the enzyme assays (panel F). SC-sterile control, n.d. - not detected.

Figure 7

Distribution of MS/MS of glucuronidated compounds among biological samples in public data across the untargeted metabolomics data from GNPS/MassIVE and Metabolights based on matches against the MS/MS resource obtained in this study. Heatmaps show the MS/MS matches and distribution of the glucuronidated standards in different tissues and biofluids with controlled vocabulary metadata available in ReDU(El Abiead et al., 2024; Jarmusch et al., 2020) in (A) rodents, positive ionization mode, (B) rodents, negative ionization mode, and (C) humans, positive ionization mode. All heatmaps are shown as the percent of samples of each tissue/biofluid with matches obtained from the repository-scale search (e.g., genistein 7-β-D-glucuronide was observed in about 25% of the lower digestive tract samples from rodents acquired in the positive ionization mode available in the repositories). The X-axes were clustered using the Bray-Curtis metric. Note that close structural isomers often exhibit similar fragmentation patterns.

Figure 8

Glucuronidated feature and GUS changes with obesity. (A) Disease distribution of MS/MS matches to glucuronidated standards among human biological samples in positive ESI public data across untargeted metabolomics data from GNPS/MassIVE and Metabolights. Heatmap shows the frequency of MS/MS matches and distribution of the glucuronidated standards in samples with controlled vocabulary disease ontology metadata (DOIDCommonName) available in ReDU and Pan-ReDU (e.g., 100% of the matches of indoxyl glucuronide to the datasets that had specified DOIDCommonName were to diabetes mellitus). (B) Volcano plot of individual GUS genes detected in HFD and NFD fed mice. Shapes of the points indicate the GUS loop class of each gene, while color indicates the statistical significance: blue=less abundant in HFD, red= more abundant in HFD. Lower case letters indicate up- or down- regulated CTD and Loop 1 GUS genes: a - 1XD42_69_MGG47451_Loop_1, b - Lachnospiraceae_unclassified_CTD, c - Unknown_unclassified_CTD, d - 1XD42_69_sp011959925_Loop_1. Significance was defined as | |Log2FC|>1.0 and p<0.05. Hypothesis testing was performed using the Kruskal-Wallis test. (C) Correlation plot of Unknown_unclassified_CTD GUS log transformed, normalized read counts and fecal aglycone m/z 220.1366 log2 transformed peak areas, with linear regression fit line (y=1.0456x+21.0364). (D) Correlation plot of Unknown_unclassified_CTD GUS log transformed, normalized read counts and urine glucuronidated feature m/z 396.1684 log2 transformed peak areas, with linear regression fit line (y=0.4358x+22.6239). Adjusted R² and p-value generated using lm in R. (E) Plot of untransformed peak areas from MSV00084112 human urine glucuronide number 8706 (m/z 396.1644), stratified by disease status. Hypothesis testing was performed using analysis of variance tests with Tukey HSD correction for multiple comparisons and statistical significance was defined at p<0.05.