Integrating Microbiome Analysis, Metabolomics, Bioinformatics, and Histopathology to Elucidate the Protective Effects of Pomegranate Juice against Benzo-alpha-pyrene-Induced Colon Pathologies

Polycyclic aromatic hydrocarbons, e.g., benzo[a]pyrene (BaP), are common dietary pollutants with potential carcinogenic activity, while polyphenols are potential chemopreventive antioxidants. Although several health benefits are attributed to polyphenol-rich pomegranate, little is known about its interaction with BaP. This study integrates histochemical, microbiomic, and metabolomic approaches to investigate the protective effects of pomegranate juice from BaP-induced pathologies. To this end, 48 Sprague–Dawley rats received, for four weeks, either pomegranate, BaP, both, or neither (n = 12 rats per group). Whereas histochemical examination of the colon indicated tissue damage marked by mucin depletion in BaP-fed animals, which was partially restored by administration of pomegranate juice, the fecal microbiome and metabolome retained their resilience, except for key changes related to pomegranate and BaP biotransformation. Meanwhile, dramatic microbiome restructuring and metabolome shift were observed as a consequence of the elapsed time (age factor). Additionally, the analysis allowed a thorough examination of fecal microbiome–metabolome associations, which delineated six microbiome clusters (marked by a differential abundance of Lactobacillaceae and Prevotellaceae, Rumincococcaceae, and Erysipelotrichaceae) and two major metabolome clusters (a sugar- and amino-acids-dominated metabotype vs. a cluster of fatty acids and hydrocarbons), with sugar alcohols maintaining a unique signature. In conclusion, using paired comparisons to minimize inter-individual animal variations allowed the dissection of temporal vs. treatment-derived variations. Microbiome–metabolome association clusters may be further exploited for metabotype prediction and gut-health biomarker discovery.


Introduction
The incidence of cancer is on the rise in developing countries [1,2] as a result of demographic shifts, climate change, pollution, and the adoption of cancer-associated lifestyle choices [3]. The global cancer burden rose to an estimated 19.3 million new cases (18.1 million cases excluding nonmelanoma skin cancer) and 10 million cancer deaths   Beanplots comparing the mucosal acidic mucin levels, expressed as the percentage of mucosal reactive mucin areas between different treatment groups (n = 6 per group). Differences between groups were analyzed for statistical significance by the Kruskal-Wallis test, followed by post-hoc pairwise Wilcoxon testing. ns = non-significant (p value > 0.05); ** = p value < 0.005.

Histopathological Alterations in Rat Colons upon Treatment
To determine relevant phenotypic alterations in the intestinal lumen and tissues, we microscopically examined hematoxylin and eosin (H&E)-stained rat colons from the four groups (Figure 1), and we quantified goblet cell-derived reactive acidic mucins by Alcian blue (pH 2.5) staining ( Figure 2).
The histopathological examination indicated that the colonic walls of both the untreated and pomegranate-fed rat groups had well-organized histological features, with intact intestinal crypts ( Figure 1A,B,E,F). The crypts had abundant goblet cells, alternated with many enterocytes, with intact nuclear and subcellular details, intact propria, and intact submucosal layers, without abnormal infiltrates. On the other hand, the colons of BaP- Beanplots comparing the mucosal acidic mucin levels, expressed as the percentage of mucosal reactive mucin areas between different treatment groups (n = 6 per group). Differences between groups were analyzed for statistical significance by the Kruskal-Wallis test, followed by post-hoc pairwise Wilcoxon testing. ns = non-significant (p value > 0.05); ** = p value < 0.005.
The histopathological examination indicated that the colonic walls of both the untreated and pomegranate-fed rat groups had well-organized histological features, with intact intestinal crypts ( Figure 1A,B,E,F). The crypts had abundant goblet cells, alternated with many enterocytes, with intact nuclear and subcellular details, intact propria, and intact submucosal layers, without abnormal infiltrates. On the other hand, the colons of BaP-fed rats had obvious histopathological changes. Focal areas of hyperplastic and dysplastic changes were common in the lining epithelium of intestinal crypts, with dark basophilic cytoplasms and oval basally situated nuclei ( Figure 1C,G). Additionally, these colons were characterized by an increased nuclear cytoplasmic ratio (N/C) with many mitotic figures, minimal differentiated mature goblet cells, and marked interstitial inflammatory cell infiltrate ( Figure 1G).
As for the rats fed with both substances (the PomBaP group), dysplastic changes in the lining epithelium of intestinal crypts were markedly reduced, with many intact enterocytes and apparent intact mature goblet cells. An inflammatory cell infiltrate was still observed in these rats, accompanied by mild congested mucosal blood vessels ( Figure 1D,H).
The highest level of reactive mucosal acidic mucin was observed in the pomegranatefed rats, followed by the untreated rats, whereas a significant reduction in reactive mucin was observed in the BaP-fed rats ( Figure 2C,E). In the PomBaP-fed rats, mucin levels were significantly higher (median~20%) than in the BaP-fed rats (post-hoc Wilcoxon test p value = 0.002), yet they were not fully restored to the mucin levels of the Pom or untreated groups (median~25%, Figure 2E).

Uneven Increase in β-Glucuronidase Activity after Four Weeks
Glucuronidases are key xenobiotic-modifying bacterial enzymes that are enriched in specific bacterial classes and have been reported to increase in the microbiota in association with some types of cancer. For this reason, we ensured that the rats were randomized for week 0 fecal glucuronidase activity, which was determined right before the feeding system was initiated ( Figure S2A indicates no significant difference between cages). After the four-week experiment, the glucuronidase activity significantly increased in the entire rat population of the experiment (paired Wilcoxon test p value = 0.0012, Figure 3A), yet there was no significant difference in week 4 glucuronidase activity between different groups ( Figure S2B). Peculiarly, the within-group increase in β-glucuronidase activity did not reach statistical significance ( Figure 3B) but tended to be significant in the Pom + BaP-fed rats (paired Wilcoxon p value 0.055), as eight out of nine rats in this group had modest to large increases in β-glucuronidase activity (paired lines between light blue dots in Figure 3A). changes were common in the lining epithelium of intestinal crypts, with dark basophilic cytoplasms and oval basally situated nuclei ( Figure 1C,G). Additionally, these colons were characterized by an increased nuclear cytoplasmic ratio (N/C) with many mitotic figures, minimal differentiated mature goblet cells, and marked interstitial inflammatory cell infiltrate ( Figure 1G).
As for the rats fed with both substances (the PomBaP group), dysplastic changes in the lining epithelium of intestinal crypts were markedly reduced, with many intact enterocytes and apparent intact mature goblet cells. An inflammatory cell infiltrate was still observed in these rats, accompanied by mild congested mucosal blood vessels ( Figure  1D,H).
The highest level of reactive mucosal acidic mucin was observed in the pomegranatefed rats, followed by the untreated rats, whereas a significant reduction in reactive mucin was observed in the BaP-fed rats ( Figure 2C,E). In the PomBaP-fed rats, mucin levels were significantly higher (median ~ 20%) than in the BaP-fed rats (post-hoc Wilcoxon test p value = 0.002), yet they were not fully restored to the mucin levels of the Pom or untreated groups (median ~ 25%, Figure 2E).

Uneven Increase in β-Glucuronidase Activity after Four Weeks
Glucuronidases are key xenobiotic-modifying bacterial enzymes that are enriched in specific bacterial classes and have been reported to increase in the microbiota in association with some types of cancer. For this reason, we ensured that the rats were randomized for week 0 fecal glucuronidase activity, which was determined right before the feeding system was initiated ( Figure S2A indicates no significant difference between cages). After the four-week experiment, the glucuronidase activity significantly increased in the entire rat population of the experiment (paired Wilcoxon test p value = 0.0012, Figure 3A), yet there was no significant difference in week 4 glucuronidase activity between different groups ( Figure S2B). Peculiarly, the within-group increase in β-glucuronidase activity did not reach statistical significance ( Figure 3B) but tended to be significant in the Pom + BaPfed rats (paired Wilcoxon p value 0.055), as eight out of nine rats in this group had modest to large increases in β-glucuronidase activity (paired lines between light blue dots in Figure 3A).

Microbiome Profiling by 16S rRNA Amplicon Sequencing
To investigate any potential effects exerted by the pro-carcinogenic BaP, anti-carcinogenic pomegranate juice, and their combination on the rat gut microbiota, we sequenced the DNA extracted from 24 samples collected from 12 cages before and after the 4-week treatment period. Each sample represented carefully collected, pooled animal feces from a different cage, and thus each treatment group (untreated, Pom, BaP, and PomBaP) included three different cages (biological replicates-with all live animals therein), while the control included week 0 samples from all 12 cages.
The raw sequencing reads were quality filtered, pre-processed, and then analyzed by MOTHUR software against the SILVA 16S rRNA database (v.123). MOTHUR analysis filtered and classified 171,078 reads, with an average of 7128 reads per sample, that were collectively assigned to 367 non-redundant taxonomic units of different taxonomic levels (including 17 phyla, 32 classes, 47 orders, 74 families, and 195 genera). Subsequently, the assigned taxa were compared among samples and sample groups at different taxonomic levels. All data were then re-filtered, statistically analyzed, and visualized by Microbiome-Analyst software [65] as follows.
Analyzed reads were normalized across samples (as a fraction of total reads) and rarefied to the minimum library size (1369), and 11,150 low-abundance features + 20 lowvariance features were removed after a low-variance filter was applied (detailed in the Materials and Methods Section).
At the phylum level, MOTHUR initially identified 17 phyla, 15 of which were retained by MicrobiomeAnalyst (as two of the initially classified phyla were filtered out for low abundance). Four phyla were the most predominant in all samples: Bacteroidetes (with mean relative abundance = 65.9% of all classified reads and 67.8% of filtered classified reads), followed by Firmicutes (with mean relative abundance = 26.4% of all classified reads and 17.4% of filtered classified reads), Tenericutes (with mean relative abundance = 5.8% of all classified reads and 14.3% of filtered classified reads), and Actinobacteria (with mean relative abundance = 1.6% of all classified reads and 0.52% of filtered ones). Eleven other phyla were found in low proportions ( Figure S3A). At the genus level, 193 genera were observed, 16 of which were the predominant ones ( Figure S3B), while 176 other genera were found in low proportions (Supplementary Data File S1).

Alterations in Fecal Microbiota Diversity across Sample Groups
The alpha diversity of the rat fecal microbiota was estimated, by different indices, at the genus level. Three out of the four calculated alpha diversity indices (observed OTUs, Shannon, and Simpson) were significantly different between control samples and each of the following treatment groups: untreated, Pom, and BaP ( Figure S4A). The same three indices were significantly different between the week 0 and week 4 samples ( Figure  S4C). However, when the significantly different control samples were excluded from the comparison, because they all belong to week 0, no other significant difference was observed between the alpha diversity metrics of any of the week 4 treatment groups ( Figure S4B).
The gut microbiota beta diversity, assessed through principal component analysis (PCA) of Bray-Curtis distances of all sequence variants (feature-level analysis), indicated no segregation between different treatment groups, except for the control group (week 0) samples ( Figure 4A). Thus, the time attribute could be the more discriminant factor. Indeed, the 24 microbial communities were separated into two distinct clusters based on age/temporal factor: week 0 vs. week 4 ( Figure 4B). At the genus level, the same clustering pattern was retained, but to a lesser extent ( Figure 4C,D).
Indeed, the 24 microbial communities were separated into two distinct clusters based on age/temporal factor: week 0 vs. week 4 ( Figure 4B). At the genus level, the same clustering pattern was retained, but to a lesser extent ( Figure 4C,D). Overall, both alpha and beta diversity analyses strongly highlighted the age or temporal factor as a major predictor of microbiome composition and diversity. Meanwhile, despite observed cage-to-cage variability ( Figure S3), the cage factor did not seem to consistently affect microbiome profiles across treatment groups, and no specific clustering by cage could be seen among samples ( Figure S5). The strong temporal factor does not discount that some changes are caused by the treatment type (i.e., group factor), which seemed less pronounced than, but distinct from, the age factor, as detailed below. Overall, both alpha and beta diversity analyses strongly highlighted the age or temporal factor as a major predictor of microbiome composition and diversity. Meanwhile, despite observed cage-to-cage variability ( Figure S3), the cage factor did not seem to consistently affect microbiome profiles across treatment groups, and no specific clustering by cage could be seen among samples ( Figure S5). The strong temporal factor does not discount that some changes are caused by the treatment type (i.e., group factor), which seemed less pronounced than, but distinct from, the age factor, as detailed below.

Alterations in Fecal Microbiome Profiles across Sample Groups
Upon comparison of the microbiome profiles among all five groups (the four treatment groups at week 4, including three samples/cages per group, as well as the control samples, including all 12 cages of week 0), 19 taxa ( Figure S6) of different taxonomic levels were deemed statistically significantly different between the groups, as calculated by linear discriminant analysis effect size (LEfSe). A shorter list of the taxa is provided in Figure 5.
Upon comparison of the microbiome profiles among all five groups (the four treatment groups at week 4, including three samples/cages per group, as well as the control samples, including all 12 cages of week 0), 19 taxa ( Figure S6) of different taxonomic levels were deemed statistically significantly different between the groups, as calculated by linear discriminant analysis effect size (LEfSe). A shorter list of the taxa is provided in Figure  5. Among these differentially abundant taxa, order Mollicutes RF9 (class Mollicutes); family Peptostreptococcaceae and genus Romboutsia (class Clostridia, order Clostridiales); and genus Veillonella were less abundant in the control group ( Figure 5A), while order Bacteroidales and family Bacteroidales S24-7 group (phylum Bacteroidetes, class Bacteroida) and family Veillonellaceae and genus Anaerovibrio were relatively more abundant in the control group ( Figure 5B).
Family Acidaminococcaceae and genus Phascolarctobacterium; Family XIII; and genus Erysipelotrichaceae UCG-003 were relatively more abundant in the PomBaP group, while genus Prevotellaceae UCG_003 was relatively less abundant in the same group ( Figure 5C). Among these differentially abundant taxa, order Mollicutes RF9 (class Mollicutes); family Peptostreptococcaceae and genus Romboutsia (class Clostridia, order Clostridiales); and genus Veillonella were less abundant in the control group ( Figure 5A), while order Bacteroidales and family Bacteroidales S24-7 group (phylum Bacteroidetes, class Bacteroida) and family Veillonellaceae and genus Anaerovibrio were relatively more abundant in the control group ( Figure 5B).
Family Acidaminococcaceae and genus Phascolarctobacterium; Family XIII; and genus Erysipelotrichaceae UCG-003 were relatively more abundant in the PomBaP group, while genus Prevotellaceae UCG_003 was relatively less abundant in the same group ( Figure 5C).

Influence of Temporal, Age-Related Factors on Fecal Microbiome Composition
Because the above analysis included two major experimental factors, namely time (i.e., week 4 vs. week 0, reflecting 28 days of lifespan) and the treatment factor (three different feeding systems, in addition to the solvent-only or untreated group), we compared all samples, lumped by sampling week.
This comparison indicated an increase in the relative abundance of some taxa with time/age. The abundance of phylum Tenericutes; classes Mollicutes and Clostridia; orders Mollicutes RF9 and Clostridiales; families Acidaminococcaceae and Peptostreptococcaceae; and genera Veillonella, Phascolarctobacterium, and Romboutsia seemed to increase at the end of week 4 ( Figure 6A). On the other hand, the abundance of phylum Bacteroidetes; class Bacteroidia; order Bacteroidales; families Prevotellaceae; Bacteroidales S24-7 group; and Veillonellaceae, genera Anaerovibrio, Prevotellaceae UCG-001, and Prevotellaceae UCG-003 seemed to decrease in week 4 samples ( Figure 6B).
Because the above analysis included two major experimental factors, namely time (i.e., week 4 vs. week 0, reflecting 28 days of lifespan) and the treatment factor (three different feeding systems, in addition to the solvent-only or untreated group), we compared all samples, lumped by sampling week.
This comparison indicated an increase in the relative abundance of some taxa with time/age. The abundance of phylum Tenericutes; classes Mollicutes and Clostridia; orders Mollicutes RF9 and Clostridiales; families Acidaminococcaceae and Peptostreptococcaceae; and genera Veillonella, Phascolarctobacterium, and Romboutsia seemed to increase at the end of week 4 ( Figure 6A). On the other hand, the abundance of phylum Bacteroidetes; class Bacteroidia; order Bacteroidales; families Prevotellaceae; Bacteroidales S24-7 group; and Veillonellaceae, genera Anaerovibrio, Prevotellaceae UCG-001, and Prevotellaceae UCG-003 seemed to decrease in week 4 samples ( Figure 6B).  Having examined the overall variability between different groups and between the week 0 and week 4 samples, we set out to resolve differences that are only contributed by the treatment factor. This required normalization to week 0 (i.e., calculating a week 4-toweek 0 abundance ratio for different microbial taxa). Because it was not mathematically possible to calculate ratios for taxa that were not detected on week 0 (i.e., having zero relative abundance values), the week 4/week 0 ratio was only calculated for taxonomic groups with non-zero values in week 0 samples (92 taxonomic units of all levels).
When the samples from the four different treatment groups of week 4 were normalized to their corresponding samples of week 0 with non-zero relative abundance, we found that family Lactobacillaceae (mainly Lactobacillus), non-Lactobacillaceae members of class Bacilli (i.e., Bacilli of unclassified order, family, or genus), and genus Allobaculum were significantly more abundant in PomBaP-fed rats, while unclassified Ruminococcacae were more abundant in the untreated rats. Additionally, unclassified genera of Veillonelleacae were significantly (ANOVA p value < 0.05) more abundant in the BaP-fed rats ( Figure 7). After repeating the same analysis with non-parametric statistical testing (Kruskal-Wallis test followed by Dunn's or paired Wilcoxon post-hoc tests), only genus Allobaculum and family Prevotellaceae (unclassified genera) were deemed statistically significantly different between groups (p value < 0.05). Allobaculum abundance was significantly higher in week 4 samples from the PomBap group and significantly lower in the week 4 Pom group than all other groups, while unclassified Prevoteallaceae were more abundant in the Pom group of week 4 than in the PomBaP group.

Alterations in Fecal Metabolome across Sample Groups
To understand the metabolic consequence of the four-week feeding scheme and its associated microbiome alterations, we further analyzed 72 fecal samples (three technical replicates of 24 samples, representing the 12 cages before and after the experiment) by gas

Alterations in Fecal Metabolome across Sample Groups
To understand the metabolic consequence of the four-week feeding scheme and its associated microbiome alterations, we further analyzed 72 fecal samples (three technical replicates of 24 samples, representing the 12 cages before and after the experiment) by gas chromatography coupled with mass spectrometry detection (GC-MS) to dissect the changes in metabolomic profiles before and after treatment, and between different treatment groups. Metabolites detected as trimethylsilylated derivatives were mostly low-molecularweight primary metabolites, such as organic acids, amino acids, fatty acids, esters, alcohols, nitrogenous compounds, phenolic acids, sterols, and sugars, in addition to aliphatic and aromatic hydrocarbons, some of which are likely derived from metabolism of pomegranate juice by gut microbiota.
One hundred and two metabolites (Table S1 and Supplementary Data File S2), out of 530 mass spectra (MS) detected peaks, were annotated from all examined samples and classified into 13 classes ( Figure 8A). The levels of the detected metabolites were normalized to that of spiked internal standard xylitol.  GC-MS analysis indicated that the rat fecal metabolome switched from a sugardominated to a fatty-acid-dominated profile ( Figures 8A and S8). On average, sugars represented the major class in samples from week 0, constituting 54% of the total detected metabolites, followed by acids (13%), whereas fatty acids represented the major class in samples from week 4, constituting 26% of the total metabolites, followed by sugars (22%). Other minor identified classes included sterols (0.44% in week 0 and 0.9% in week 4), esters (0.47% in week 0 and 0.64% in week 4), and aliphatic hydrocarbons (0.5% in week 0 and 1.15% in week 4 (Table S2)).
As for the major changes in the levels of metabolite classes each of week 4 treatment groups compared to the untreated group two classes seemed to diminish in all groups: (i) acids diminished by 0.29, 0.26, and 0.31-fold in the Pom, BaP, and PomBaP groups, respectively; and (ii) aromatic hydrocarbons diminished by 0.42, 0.47, and 0.53-fold in the Pom, BaP, and PomBaP groups, respectively.
The following classes were enriched in all treatment groups, compared to the untreated group (Table S1): (i) alcohols increased by 1.7, 6.6, and 1.4-fold in the Pom, BaP, and PomBaP groups, respectively; (ii) amino acids increased by 2.2, 1.3, and 2.8-fold in the Pom, BaP, and PomBaP groups, respectively; (iii) phenolic acids increased by 2.5, 2.1, and 1.8-fold in the Pom, BaP, and PomBaP groups, respectively; (iv) sugars increased by 1.6, 1.4, and 1.1-fold in the Pom, BaP, and PomBaP groups, respectively; and (v) sugar alcohols increased by 3.6, 2.3, and 8.8-fold in the Pom, BaP, and PomBaP groups, respectively. Fatty acids were enriched in the Pom and BaP groups (by 1.15 and 1.16-fold, respectively) but diminished in the PomBaP one (by 0.5-fold).

Unsupervised and Supervised Multivariate Data Analysis of GC-MS Datasets
To further classify the treatment groups, we used PCA and orthogonal projection leastsquares discriminant analysis (OPLS-DA) to determine the impact of different treatments on rat metabolome profiles in an untargeted manner.

Clustering of the Entire Dataset into Two Separate Clusters by PCA and OPLS-DA
We started by applying PCA to the entire dataset (all samples and all metabolites included) to investigate potential clustering patterns and determine possible metabolome heterogeneity among all examined samples. Nine components accounted for 90.4% of the total variance (R 2 ), with the first two components (PC1/PC2) accounting for 65% of the variance (score plot in Figure 8C). All samples could be segregated into two distinct clusters, distributed along PC1. As observed with the microbiome data, these two major clusters were clearly set apart based upon the age/temporal factor (week 0 vs. week 4, Figure 8B-D).
To determine the key metabolites responsible for the sample segregation, we applied loading plots to the PCA models. Among detected metabolites, sugars (glucose, arabinose, and glycerol), fatty acid esters (1-monooleoylglycerol), and aromatic hydrocarbons (p-xylene) accounted for most of the observed variability along PC1 ( Figure 8C). Metabolites with extreme positive p [1] values were more abundant in week 4 samples and contributed to their separation, whereas metabolites with extreme negative p [1] values were more abundant in week 0 samples and led to their separation to the left of the graph.
Considering that sugars accounted the most for the observed separation between week 0 and week 4 samples ( Figure 8C), we repeated the PCA analysis after the exclusion of all sugars. The same segregation pattern was largely retained in this analysis, in which PC1/PC2 explained 55% of the variance ( Figure S9). However, while week 0 samples were clearly clustered, week 4 samples were intriguingly divided into a major cluster, which included most samples, and a minor one, which was enriched in samples from the Pom-treated group ( Figure S9A). Such a pattern was confirmed by hierarchical clustering analysis, which separated the non-sugar metabolome profiles into three clades, including a week 4 Pom-enriched clade ( Figure S9C).
After the exclusion of sugars, the corresponding loading plots uncovered another set of metabolites that contributed to the variability along PC1 and, hence, segregation. This set of metabolites included acids (hydroxybutyric acid, hydroxybenzoic acid), the amino acids valine, the alcohol 1,4-butanediol, and the previously detected fatty acid esters (1-monooleoylglycerol) and aromatic hydrocarbons (p-xylene, Figure S9B).
In agreement with the PCA results, a supervised OPLS-DA model, applied to the entire dataset and explaining 48% of the total variance, confirmed the age-based segregation of the samples ( Figure S10A). Likewise, as sugars constituted the major variable separating week 0 from week 4 samples ( Figure 8C), we re-applied the same model after the exclusion of sugars to uncover other possible contributors to the observed age/time-based segregation pattern ( Figure S10C).
For OPLS-DA models, an S-loading plot was used to visualize the covariance (p) and the correlation (pcor) between the variability-generating metabolites. Like with PCA, sugars (glucose, arabinose, and glycerol), the fatty acid ester 1-monooleoylglycerol, and the aromatic hydrocarbon p-xylene accounted for most of the observed variability along PC1 ( Figure S10B). Yet, unlike PCA, hydroxybutyric acid, the fatty acid myristic acid, and the sugar alcohol maltitol added to observed variability along PC1 ( Figure S10B).
After the exclusion of sugars from the OPLS-DA model, like with PCA, 1-monooleoylglycerol, hydroxybutyric acid, p-xylene, 1,4-butanediol, and the amino acid valine accounted for most of the observed variability along PC1. Yet, unlike PCA, fatty acids and fatty acid esters (1-monopalmitin and myristic acid), 3-hydroxyphenylpropionic acid, and the amino acid L-leucine added to the observed variability along PC1 ( Figure S10D).
Of note, because of the complexity of data and multiple compared variables, we conducted week 0 vs. week 4 analysis in subsets based on the assigned treatment groups (Figures S11 and S12). An OPLS-DA model and S-Plot that were applied to each of week 0-week 4 pairs singled out sugars as the dominant class of discriminatory metabolites ( Figure S11). After the exclusion of sugars, several metabolites seemed to contribute to the distinction of week 0 samples to include succinic acid (in the untreated and Pom groups); linoleic acid (in the BaP and PomBaP groups); isoleucine (enriched in the Pom group); and uracil, palmitic acid, and 9-octadecenoic acid -2-propyl ester (enriched in the PomBaP group). Meanwhile, lactic acid was associated with week 4 samples of the PomBaP group ( Figure S12).

Lack of Segregation between Week 0 Sample Groups
PCA was also applied to the normalized metabolome profiles from week 0 samples to investigate any metabolomic heterogeneity among these samples from control rats before any treatment was applied. As expected, the entire pool was not segregated into any subclusters ( Figure S13). A loading plot for this analysis confirmed the absence of major discriminatory metabolites among week 0 samples, except for some individual metabolites marking the natural biological variability among rats. These metabolites included lactic acid, as well as glucose/glucopyranose and p-xylene, which were associated with a couple of outliers ( Figure S13B).

Alterations in Fecal Metabolomes in Response to Daily Ingestion of Pom, BaP, and PomBaP
Like in analyzing microbiome data, it was important to determine treatment-specific factors that are independent from age/temporal changes. Thus, metabolite level values from different treatment groups of week 4 were normalized to their corresponding values in week 0 samples. When PCA, as an unsupervised model, was applied to this double-normalized dataset of week 4 samples, 77.5% of the variance was explained by five components, with PC1/PC2 accounting for 48.6% of the total variance. The samples could be separated into two clusters distributed along PC1, a cluster representing the untreated group and another cluster representing the rest of the treatment groups. (Figure 9A). Sugar alcohols (maltitol, myoinositol, pinitol), the amino acid alanine, and propionic acid accounted for this separation ( Figure 9B).

Alterations in Fecal Metabolomes in Response to Daily Ingestion of Pom, BaP, and PomBaP
Like in analyzing microbiome data, it was important to determine treatment-specific factors that are independent from age/temporal changes. Thus, metabolite level values from different treatment groups of week 4 were normalized to their corresponding values in week 0 samples. When PCA, as an unsupervised model, was applied to this doublenormalized dataset of week 4 samples, 77.5% of the variance was explained by five components, with PC1/PC2 accounting for 48.6% of the total variance. The samples could be separated into two clusters distributed along PC1, a cluster representing the untreated group and another cluster representing the rest of the treatment groups. (Figure 9A). Sugar alcohols (maltitol, myoinositol, pinitol), the amino acid alanine, and propionic acid accounted for this separation ( Figure 9B).  As a supervised model, OPLS-DA accounted for 59.1% of the total variance, and the examined dataset could also be separated into three clusters along PC1, which were more distinct than in the PCA, as expected from supervised modeling as such. The clusters included an untreated group cluster, a PomBaP group cluster, and a cluster combining the Pom and BaP groups ( Figure 9C). The same metabolites (maltitol, myoinositol, pinitol, alanine, and propionic acid) accounted for this separation, in addition to leucine and octadecan-1-ol ( Figure 9D).
To unveil other metabolomic differences between week 4 treatment groups, we implemented a set of OPLS-DA models for pairwise comparison between these groups ( Figure 10). The model comparing the untreated vs. Pom groups (R 2 = 98.8, Figure 10A) indicated propionic acid, linoleic acid, maltitol, and alanine as the principal contributors to the observed variability ( Figure 10G). The model comparing the untreated vs. BaP groups (R 2 = 98.6%, Figure 10B) suggested propionic acid, maltitol, and ethylene glycol as the major variables ( Figure 10H). As for the untreated vs. PomBaP model (R 2 = 96.5%, Figure 10C), it suggested propionic acid, maltitol, myoinositol, and pinitol as the most variable metabolites along PC1 ( Figure 10I).
In these models, 1-monooleoylglycerol, 9-octadecenoic acid -2-propyl ester, 4hydroxyhydrocinnamic acid, 5-hydroxyindoleacetic acid, oleic acid, and adenosine accounted for most of the variability between the Pom and PomBaP groups ( Figure 10J); 1-monooleoylglycerol and ethylene glycol accounted for most of the variability between the BaP and PomBaP groups ( Figure 10K); and adenosine, ethylene glycol, benzenepropanoic acid, and 4-hydroxy-3-methoxy were the most variable metabolites along S-Plot PC1 between the Pom and BaP groups ( Figure 10L).

Microbiome-Metabolome Correlations
Having separately analyzed compositional and metabolic variations in the fecal microbiota, we investigated the potential metabolome-microbiome profile correlation(s) of all 24 samples, regardless of sampling time or treatment group (Figures 11 and S14). In other terms, we determined which metabolites are more frequently associated with which microbial taxa in a statistically significant manner (with a Spearman correlation coefficient, r s , cutoff of ±0.5, p value < 0.05, and a taxon's relative abundance threshold of 0.05% in at least one sample-Supplementary Data File S3). Overall, we found several amino acids to be more often correlated with Bacteroidetes (especially members of Bacteroidales S24-7 family, with r s = 0.783), while amino sugars were correlated with class Bacilli of the Firmicutes phylum (mainly represented by genus Lactobacillus, r s = 0.547) and genus Anaerovibrio of class Negativicutes (r s = 0.583). On the other hand, sugar acids were negatively correlated with phylum Firmicutes (mostly contributed by genus Veillonella and class Clostridia, with r s = −0.557 and −0.715, respectively ( Figures 11A,B and S14)). As for the major class 'sugars', they were positively correlated with Bacteroidetes (r s = 0.643) in general, especially class Bacteroidia (r s = 0.742) and family Bacteroidales S24-7 (r s = 0.814), in addition to some members of phylum Firmicutes, especially Anaerovibrio (r s = 0.730) and Lactobacillus (r s = 0.543); conversely, sugars were negatively correlated with Family XIII (r s = −0.773), and genera Anaerovorax and Ruminococcaceae UCG-014 (r s = −0.808 and −0.807, respectively). cated propionic acid, linoleic acid, maltitol, and alanine as the principal contributors to the observed variability ( Figure 10G). The model comparing the untreated vs. BaP groups (R 2 = 98.6%, Figure 10B) suggested propionic acid, maltitol, and ethylene glycol as the major variables ( Figure 10H). As for the untreated vs. PomBaP model (R 2 = 96.5%, Figure  10C), it suggested propionic acid, maltitol, myoinositol, and pinitol as the most variable metabolites along PC1 ( Figure 10I).    Dominant members of phylum Firmicutes were divided into two contrasting patterns in how they correlate with metabolite levels: Most members of class Bacilli, e.g., family Lactobacillaceae and genus Lactobacillus but not family Streptococcaceae, had a positive correlation with amino acids, sugars, alcohols, and sugar acids (rs = 0.428, 0.543, 0.525, and 0.622, respectively), but a negative correlation with fatty acids and aliphatic hydrocarbons (rs = −0.677 and −0.394, respectively); by contrast, class Clostridia, including families Ruminococcaceae, Peptostreptococcaceae, and Family XIII, had metabolite correlations that are more similar to class Mollicutes of the phylum Tenericutes, which tended to correlate with sterols, fatty acids and fatty acid esters, and aliphatic hydrocarbons ( Figure 11). Dominant members of phylum Firmicutes were divided into two contrasting patterns in how they correlate with metabolite levels: Most members of class Bacilli, e.g., family Lactobacillaceae and genus Lactobacillus but not family Streptococcaceae, had a positive correlation with amino acids, sugars, alcohols, and sugar acids (r s = 0.428, 0.543, 0.525, and 0.622, respectively), but a negative correlation with fatty acids and aliphatic hydrocarbons (r s = −0.677 and −0.394, respectively); by contrast, class Clostridia, including families Ruminococcaceae, Peptostreptococcaceae, and Family XIII, had metabolite correlations that are more similar to class Mollicutes of the phylum Tenericutes, which tended to correlate with sterols, fatty acids and fatty acid esters, and aliphatic hydrocarbons ( Figure 11). Some metabolites were specifically correlated with a small number of taxa, e.g., octadecenoic acid, its ester, and eicosanoic acid were correlated (r s = 0.616, 0.612, 0.623, respectively) with the rare phylum Fusobacteria. Lactic acid, expectedly, was highly correlated with phylum Firmicutes, which includes the well-known lactic acid bacteria. Glycolic acid and glycine clustered and were positively correlated with Verrucomicrobia ( Figure S14), while pinitol, which also partly clustered with glycolic acid and glycine, was positively correlated with Verrucomicrobia, Erysipelotrichia, and Firmicutes in general (r s = 0.44, 0.608, and 0.481, respectively). Fatty acids and fatty acid esters, especially the propyl esters of stearic acid and 9-Octadecenoic acid, were positively correlated with phylum Tenericutes and its class Mollicutes, as well as Proteobacteria and Fusobacteria, but were interestingly most correlated with unclassified bacteria (Figures 11 and S14). A full list of correlation coefficients is provided in Supplementary Data File S3.
In addition to metabolites determined by GC-MS, the measured fecal β-glucuronidase activity was positively correlated with the relative abundance of class Clostridia, especially Family XIII, genera Veillonella and Holdmanella, and family Ruminococcaceae (with correlation coefficients r s = 0.56, 0.62, 0.65, 0.66, and 0.57, respectively) and negatively correlated with Bacteroidetes (r s = −0.56) and genus Anaerovibrio of the Negativicutes (r s = −0.53).
Of interest, when we examined the correlation among bacterial taxa based on their co-abundance with fecal metabolites (Figures 11C and S15), two major clusters of bacterial taxa were distinguishable ( Figure 11C): one cluster includes Bacilli, Bacteroidetes, and Negativicutes, while the other includes Clostridia and most Firmicutes, Verrucomicrobia, and Tenericutes. At the family and genus levels, these two clusters could be further divided into six ( Figure 12 and Table 1). Table 1. Microbiome clusters (visualized in Figure 12).

Cluster
Major Members Likewise, metabolic classes were split into two major clusters: one included amino acids, sugars, amino sugars, and alcohols, while the other included fatty acids and their esters, as well as aromatic and aliphatic hydrocarbons. Sugar alcohols stood out uniquely, with correlation patterns that spanned the two clusters ( Figure 11D,E). Correlation patterns of individual metabolites are provided in the Supplementary Materials ( Figure S16).

Discussion
Humans and animals are continually and increasingly exposed to environmental toxicants, many of which are potential carcinogens. Among these toxicants are polycyclic aromatic hydrocarbons (PAHs), which are common air pollutants and are enriched in grilled, fried, and processed foods [21][22][23]. A prototypic representative of these toxicants is benzo[a]pyrene (BaP) [24,25], which is one of the two main chemicals studied in this work.
As a reaction to the alarming levels of pollutants, especially carcinogens, in food and the environment, several calls have been made to return 'back to nature' by seeking alternatives to procarcinogenic foods and toxicants or at least protection from their harmful effects. Among the most popular natural products are polyphenols, and their representative plant source in this study is pomegranate juice.
Plenty of polyphenol-rich foods have been examined for anticarcinogenic properties against BaP-induced colon pathologies in vivo. These food components include resveratrol [66], olive oil [67], silver Moringa oleifera leaves [68], and curcumin [69]. However, little experimental evidence is available in the literature on the role of pomegranate as a chemoprotective agent in the presence of BaP.
Great attention has lately been paid to the role of the human microbiota, the assemblage of microorganisms residing in different living and non-living environments, and among the key described roles of the microbiota is its role in modulating xenobiotics (toxicomicrobiomics), sometimes increasing their toxic effects and sometimes detoxifying them [70,71].
While experimentally investigating the causality of human-microbiome interactions is usually neither ethical nor feasible, animal models are powerful tools for establishing causality, even if not fully representative of the human host. They allow the use of proper controls and allow longitudinal studies to be performed in a relatively short time, without the obstacles and limitations associated with human studies.
Here, we used a well-controlled rat model to investigate the changes in the intestinal tissue, microbiome, and metabolome over time, in response to the ingestion of two chemicals with procarcinogenic and anticarcinogenic potential.
The major phenotypic changes that we observed in the experimental model are histopathological, while animal weights remained unchanged. In addition, we measured the glucuronidase activity in feces as a marker for xenobiotic-modifying potential, and we ensured that rats were randomized for this activity, which significantly increased after the four-week experiment, but that increase was similar across all treatment groups, and no significant increase was observed in a particular group.

Histopathological Alterations in Rat Colons
Our histopathological examination indicated significant observable inflammation and damage to the colonic tissues of rats after a daily ingestion of BaP for 28 days, which was not detected at all in Pom-fed rats and was partly restored by the addition of Pom to BaP (Figures 1 and 2). The damage was manifested by signs of inflammation, mitosis, dysplasia, and a significant decrease in reactive mucin.
Mucins are produced and secreted by goblet cells and represent a source of nutrition to beneficial commensals while protecting the colon epithelium from pathogenic bacteria by forming a covering mucus layer [72]. Pomegranate's positive impact on mucin is well documented. For example, supplementation of the basal diet administered to male Sprague-Dawley rats with 10% pomegranate powder significantly increased the mucin content in tissues that had aspirin-induced ulcers by 52.7%, compared to untreated animals [73]. Pomegranate extract significantly decreased the number of mucin-depleted foci (MDF) in the colons of rats fed a carcinogenic dark cooked meat with nitrite, oxidized (DCNO)-based diet [74].
When a pomegranate mesocarp decoction was administered, at a dose of 50 mg/kg/day for 6 weeks, to a rat model that spontaneously developed colon adenomas, it significantly reduced the number of MDF [75]. In a study on a chronic trinitrobenzene sulfonic acid (TNBS)-induced colitis Wistar rat model, dietary supplementation of pomegranate extract to the rats for four weeks attenuated the morphological signs of cell damage in the colonic mucosa and reversed mucin depletion in the Alcian-blue-positive goblet cells [76]. In another TNBS-induced colitis model in male Wistar rats, ingestion of ellagic acid (10-20 mg/kg)-the main constituent of pomegranate-by oral gavage reduced the inflammatory cell count and morphological signs of cell damage. It also caused the colon epithelium to remain intact and increased the amount of mucin stained by Alcian blue in colon mucosa [77].

Temporal/Age Effects
The most striking finding from microbiome and metabolome analyses is that week 4 vs. week 0 differences are the major factor influencing both microbial taxonomic composition and metabolome profiles. Given that the samples were collected in a similar way, aliquoted, and immediately frozen, and that all DNA extraction, amplicon sequencing, and metabolome profiling by GC-MS were performed simultaneously and with all samples randomized at every analysis step, we believe that this strong temporal factor reflects real biological effects rather than technical variability.
This factor could be due to natural aging or could simply reflect stochastic, temporal variation of the microbiota [78]. The rats were already out of adolescence [79] and were thus neither in a transitional developmental stage (as confirmed by their weights that did not significantly increase during the four weeks of the experiment, Figure S1). They cannot be described to be aged either, because the experiment ended when the rats were 12-13 weeks old [79,80]; however, four weeks in an adult rat's lifespan are approximately equivalent to 2.5 years in human life [80,81].
Importantly, although rats can live up to 3.5 years [80], we opted to use young adult rats (7-12 weeks from the time they were obtained until the end of the experiment) to rule out any increased cancer incidence due to the age factor [82].
Both age and temporal variations have been previously described to significantly affect microbiome composition in humans and other mammals [83,84]. Additionally, there is no specific reason to assume that these differences are related to dietary or laboratory environmental factors because the rats were acclimated to the laboratory environment and to a standard diet that was unchanged throughout the experimental period.
Among the key observed age/temporal differences is a general increase in alpha diversity (both richness and evenness, Figure S4), reflected in an increase in the Shannon diversity index as well, which takes both diversity richness and evenness into account ( Figure S4). This diversity increase was curiously accompanied by two opposing observations: one is the decrease in relative abundance of three major taxa ( Figure 6B): (i) Bacteroidetes, with the two major families of Bacteroidales S24-7 group and Prevotellaceae, (ii) Veillonellaceae, and (iii) Anaerovibrio; the other is the relative increase in unclassified bacteria (13.3% week 0 vs. 19.8% week 4, paired Wilcoxon test p value < 0.05). However, the drop in relative abundance of the three major groups mentioned above was counterbalanced by an increase in the relative abundance of a dozen other taxa, including Tenericutes, Clostridia, and Acidaminococcaceae ( Figure 6A). Family Veillonellaceae was split in terms of its temporal variations: The abundance of genus Veillonella increased ( Figure 6A), while that of genus Anaerovibrio decreased, influencing the entire family's relative abundance ( Figure 6B).
The impact of age differences, even if small, is not unprecedented. Zhang et al. [85] found significant aging-associated changes in gut microbiota profiles and serum metabolites: Phylum Firmicutes was positively correlated, while Bacteroidetes, Proteobacteria, and Actinobacteria were negatively correlated with age [85]. In another study, age was reported as the major factor, significantly altering the fecal microbiota composition of Zucker rats, outweighing genetic factors [86]. The Firmicutes-to-Bacteroidetes ratio and the relative abundance of families Bacteroidaceae and Peptostreptococcaceae decreased, while the relative abundance of families Ruminococcaceae and Bifidobacteriaceae increased with age. The same study reported the cage environment as a key factor affecting microbiome variation [86]. The differences between some reported age effects (e.g., Peptostreptococcaceae) and those that we report should not be alarming as some taxa may be naturally fluctuating, which is not necessarily related to biological aging.
Another relevant example of age effects is that antibiotic-induced microbiota depletion prior to repeat mild traumatic brain injury (RmTBI) reportedly led to chronic changes in microbiome composition and diversity, as well as metabolome profile, in both adolescent and adult Sprague-Dawley rats, with adolescents exhibiting greater changes [87].
Animal age also has a major impact on metabolomic profiles. Both multivariate data analyses (PCA and OPLS-DA) indicated a general decrease in sugar levels (mainly glucose, arabinose, and glycerol) in week 4 samples. Valine, L-leucine, and 1,4 butanediol levels also decreased in week 4 samples, and their delineation was made possible when all sugars, the dominant metabolites, were excluded from analysis models.
Reciprocally, p-xylene and the fatty acid 1-monooleoylglycerol levels increased in week 4 samples as detected by all models (PCA and OPLS-DA, sugars included and excluded). However, the exclusion of sugars allowed the detection of more changes, as short chain fatty acids (SCFAs) (i.e., hydroxybutyric acid, aromatic acids i.e., hydroxybenzoic acid, 3-hydroxyphenylpropionic acid) and fatty acids (i.e., myristic acid and 1-monopalmitin) increased in week 4 samples. Excluding a dominant factor from the analysis is an efficient technical maneuver to overcome its masking effect, especially in similar cases with big datasets and multiple time points/analysis groups.
Song et al. [88] found a significant difference in fecal fatty acid levels between CRC patients and healthy controls. However, the authors also noted a significant difference between the mean age of the two groups, which suggests that fatty acid metabolome might have been affected by age [88].
In general, we resorted to using PCA and OPLS-DA as dimensionality reduction analysis tools because of the complexity of the acquired data (72 chromatograms of 24 samples in triplicates, with 530 detected metabolites, including 102 that were fully identified and classified).
Another way to sort out the signal from noise, to resolve dominant factors (e.g., age here), and to detect smaller true-positive changes in these large datasets is subset analysis. When we conducted week 0 vs. week 4 analysis in subsets based on the assigned treatment groups, we detected an increase in 1-monooleoylglycerol with all assigned treatments (i.e., it is an age-dependent change and has no direct relation to the type of food/treatment ingested by the rats, Figure S12), while p-xylene only increased in the week 4 Pom and week 4 PomBaP groups. The increase in 3-hydroxyphenylpropionic acid level was observed in the week 4 Pom and week4 BaP groups, whereas the increase in myristic acid was observed in the week 4 untreated and week 4 PomBaP groups. In the literature, p-xylene was found among the detected phenolic metabolites in the feces of the pomegranate-peel-extracttreated group of beef calves [89].
As for the major changes in metabolite levels in the week 4 treatment groups, another technical maneuver was to apply a set of OPLS-DA models for pairwise comparison between the double-normalized datasets of the week 4 groups. Compared to the untreated group, all groups had higher levels of maltitol, the Pom group had higher levels of linoleic acid and alanine, while the BaP group had higher levels of ethylene glycol ( Figure 10). Propionic acid appeared to diminish in all groups in comparison to the untreated group, an apparent decrease that possibly arose inadvertently as the propionic acid level was unusually lower in week 0 animals that were assigned to the untreated group than other animals, which led to a higher week 4-to-week 0 ratio (Supplementary Data File S1).
Comparing the BaP or Pom groups head-to-head with the PomBaP (double-fed) group revealed interesting differentially abundant metabolites: Oleic acid, adenosine, 4hydroxyhydrocinnamic acid, 5-hydroxyindoleacetic acid, and 9-octadecenoic acid-2-propyl ester had higher levels in the Pom group than the PomBaP group, while ethylene glycol was enriched in the BaP group, in comparison with PomBap. Ethylene glycol was also enriched in the BaP group compared to the Pom group ( Figure 10).
Both 4-hydroxyhydrocinnamic acid and 3-hydroxyphenylpropionic acid are lowmolecular-weight phenolic acids that have been described to result from the gut microbiota metabolism of polyphenols [90,91]; thus, their association with the Pom group is reasonable and presents potential markers for pomegranate consumption in stool samples. Similarly, 5-hydroxyindoleacetic acid, which is a tryptophan metabolite, was found to significantly increase after dietary supplementation of magnolol polyphenols to a dextran sulfate sodium (DSS)-induced ulcerative colitis mouse model [92], and supplementation of the basal diet of male Ross-308 chicks with a mixture of pomegranate and onion extracts significantly increased oleic acid level [93].

Treatment Effects
Despite the quite strong signal contributed by the temporal factor, the paired design of the study allowed measuring changes that are associated with the experimental treatment with pomegranate, BaP, or both. For proper minimization of the temporal effect, a week 4-to-week 0 abundance ratio was calculated for different microbial taxa in every week 4 sample. This normalization was only technically possible for samples that had nonzero abundance value in every week 0 cage because division by zero is not possible. However, this data filtration step had a minor effect on the results because taxa with zero abundance values in week 0 samples had limited abundance in week 4 samples, in most cases (Supplementary Data File S1).
As for the major differences between treatment groups of week 4, unclassified genera of Prevotellaceae were significantly more abundant in the Pom-fed rats, while Allobaculum (family Erysipelotrichaceae) was less abundant in the same group. Unclassified genera of Veillonelleacae were significantly more abundant in the BaP-fed rats, while Lactobacillaceae (mainly Lactobacillus), non-Lactobacillaceae members of class Bacilli (Bacilli of unclassified order, family, or genus), and genus Allobaculum were significantly more abundant in PomBaP-fed rats.
Family Erysipelotrichaceae (including genus Allobaculum) has been associated with inflammatory gastrointestinal disorders and higher levels of proinflammatory tumor necrosis factor [94,95]. Compared to the healthy controls, an increased abundance of Erysipelotrichaceae was observed in colorectal cancer patients [96], and an increased abundance of genus UCG-003 of family Erysipelotrichaceae was observed in patients with inflammatory bowel disease [97]. A significantly higher level of Erysipelotrichaceae was found in the tumor group of an animal model of 1,2-dimethylhydrazine-induced colon cancer [98]. Urolithin A (a gut-microbiota-mediated metabolite of pomegranate) could significantly decrease the enriched Erysipelotrichaceae in mice in which intestinal damage was induced with ionizing radiation [99].
An unusually high abundance of family Veillonelleacae was observed in patients with colon cancer [100,101] and Crohn's disease [102]. Genus Veillonella was reportedly more abundant in patients with Hirschsprung's-disease-associated enterocolitis [103], colon cancer [104], and ulcerative colitis and primary sclerosing cholangitis with concomitant ulcerative colitis [105] than in healthy controls (in all cases). The abundance of family Prevotellaceae (mainly Prevotella) significantly increased after four weeks of pomegranate extract ingestion by healthy volunteers [106]. Moreover, in an experimental autoimmune encephalomyelitis mouse model, family Prevotellaceae was enriched in pomegranate-extractfed mice, compared to the control mice [107].
All the previous examples have major similarities with microbiome results reported here; however, to the best of our knowledge, no studies have reported microbiome/metabolome changes caused by a combination of BaP with pomegranate.

Microbiome-Metabolome Correlations
In general, microbiome composition per se is not as biologically meaningful as metabolome variations because it is the metabolic activity of microorganisms that more profoundly affects the host biology-an observation that was highlighted over a century ago [108]-and because several bacteria from distant taxonomic groups may have similar metabolism, and vice versa [64]. However, attempts to integrate microbiome and metabolome variations across the samples provided strong correlations between subgroups of bacteria and subgroups of metabolites ( Figures 11A,B and S14).
Even more interestingly, these associations allowed better understanding of internal correlations within metabolite classes and microbial taxa. Two distinct clusters of bacterial taxa were identified, and these were further divided into six microbiome clusters (Table 1, Figures 11C and S15), which were distinct in their co-abundance with metabolic classes ( Figure 11E). Similarly, metabolites were split into 'metabotypes': a small cluster of fatty acids, as well as aliphatic and aromatic hydrocarbons vs. a larger cluster including sugars, amino acids, alcohols, and nitrogenous compounds. Interestingly, sugar alcohols were the most distinct in their microbial co-abundance patterns (Figures 11D and S16).
In the literature, fecal amino acids, including glycine, were found to be positively correlated with families Lactobacillaceae (mainly Lactobacillus) and Verrucomicrobiaceae in Wistar rats [109]. Mao et al. [110] found that amino acids isobutyrate and isovalerate were positively correlated with Anaerovibrio in fecal samples collected from cows [110]. In a human study of cirrhotic patients, fatty acids (hexanoic acid and heptanoic acid) and the fatty alcohol hexadecanol were positively correlated with Tenericutes [111].
It is important not to interpret the microbiome-metabolome co-occurrence and coabundance patterns, which we and others reported, as a direct result of microbial metabolic processes. The main reason for this precaution is that fecal metabolomics measures both host-derived and microbial metabolites. Additionally, it is impossible to sort out cause and effect in these cases, i.e., whether the abundance (or scarcity) of some metabolites allowed specific bacterial types to flourish or the abundance of certain microbes led to the production (or consumption) of specific metabolites. Taking microbial interactions into account even further complicates the interpretation of these associations. For example, some microbes may abundantly produce metabolites (e.g., fucose, rhamnose, or lactic acid) that are immediately consumed by other species, so the net outcome will seem as if those metabolites are diminished in the presence of the bacteria that is expected to produce them. Despite these analysis challenges, microbiome-metabolome positive and negative correlation data are extremely useful, and with the accumulation of such data from various studies, models for microbiome/metabolome biomarkers are expected to abound in the near future, leading to the proposal of rapid assays for predicting different health conditions affected by the microbiota.

Other Phenotypic Alterations
In addition to the major dependent variables measured in this multivariate study (i.e., histopathological changes, mucin levels, microbiome composition and diversity metrics, and metabolite levels), animal weight and β-glucuronidase activity were measured and used for randomization to minimize their confounding effect on experimental factors.
Overall, we observed no significant change in weight in any experimental group, or between the start and end of the experiment. Other studies reported similar results with either BaP or pomegranate. For example, no significant weight change was reported between the control and BaP-fed (dose = 150 µg/kg/day) rats after 30 days of ingestion [112]. The weights of Sprague-Dawley rats did not significantly differ between the control group and either of the low or high-dose pomegranate-peel-extract-fed rats after 12 weeks of extract ingestion [113].
β-glucuronidase is one of the most studied families of microbially derived enzymes involved in the breakdown and metabolism of glycoconjugates that reach the colon intact. It deconjugates hepatically glucuronidated metabolites, reversing phase II metabolism carried out on endo-and xenobiotics and reactivating their aglycone form (which can be excreted), undergoing further enterohepatic recirculation or further degradation and biotransformation by the gut microbiota (in case of polyphenols [114]), thus regulating the pharmacokinetics as well as the level of metabolites of both polyphenols (pomegranate) and xenobiotics (pro-carcinogenic BaP), and affecting their potential favorable and/or adverse effects on host health [115][116][117].
In general, β-glucuronidase activity was described in the four major gut-associated bacterial phyla: Bacteroidetes, Firmicutes, Proteobacteria, and Actinobacteria [118]; however, within each phylum, there is a wide pattern of presence/absence of these enzymes among families and genera. Haiser and Turnbaugh [119] highlighted major differences among Firmicutes and Bacteroidetes, based on genomic analysis [119]. In human feces, high β-glucuronidase activity was observed in isolates belonging to Clostridium spp. [120] and in Roseburia and Faecalibacterium (both belonging to class clostridia) [121]. Clostridia were among the positively correlated taxa with β-glucuronidase activity in this work, as detailed in the Results Section.
A general, but non-group-specific, increase in fecal β-glucuronidase activity between week 0 and week 4 samples was observed in this study. Like with microbiome/metabolome age-related changes, it is hard to determine whether the changes are caused by the aging process or are stochastically affected by transient microbial restructuring. Reports in the literature are discrepant. For example, a comparison between adolescent and adult groups of C57BL/6 mice found no significant difference in enzyme activity in males but a significant reduction in activity in females [122]. Another study in a colon cancer BALB/c mouse model reported an increase in β-glucuronidase level with age in the control group [123]. In male Fischer rats, β-glucuronidase enzyme activity increased with age during the juvenile development stage [124] and in both the grain-based and meat-based diet groups [125], but it did not significantly change between three age groups: young (2-4 months), middle-aged (12-14 months), and old (22-24 months) [126].

Limitations of This Study
A study with such complexity and multiple interdependent variables can always be improved by a larger sample number and deeper or more comprehensive shotgun sequencing; however, for the purpose of defining key taxa and key metabolites while minimizing individual variations via the paired design (week 0 and week 4 sampling), this study achieved its goals. For 16S rRNA amplicon analysis, the depth of iSeq output is adequate for microbiome studies; however, the amplicon length (150-nucleotide read length, i.e., 300 nucleotides for a read pair) limits taxonomic resolution to the genus level, and no claims can be made about species distribution (even if the analysis pipelines may assign some species). Future studies should take advantage of the surging long-read sequencing technologies to have a higher resolution taxonomic assignment. In addition, we chose GC-MS to cover a wide range of microbial metabolites and compare across samples; however, detection of high-molecular-weight phenolics could not be achieved by GC-MS, for example. Future work will resort to both GC-MS and liquid chromatography mass spectrometry (LC-MS) for better profiling of biotransformed secondary metabolites of pomegranate.
An obvious limitation of any animal model is that laboratory animals do not exactly mimic the human system. In addition, the used Sprague-Dawley rats are neither isogenic to minimize the genetic variability nor wild to be diverse enough in their microbiota composition. However, starting with such experiments is pivotal to establishing working hypotheses that can be further validated in various genetic backgrounds, different animals, and eventually in humans.
In studies with this number of variables and groups, it is quite important to distinguish between statistical and biological significance. The number of samples per group and the presence of multiple treatment groups (untreated, Pom, BaP, PomBaP) are certainly technical hurdles against strong statistical significance; however, using this design is key for ruling out spurious, false-positive findings. For example, had we compared the different treatment groups at the end of the experiment without normalizing to week 0 values (to control for inter-individual and inter-cage variability), we would have ended up with a completely different set of results.
Meanwhile, it is important to note that many true differences are subtle and can be biologically significant, even in observed numbers or levels that did not reach statistical significance by hypothesis-testing methods. A clear example is the appearance/disappearance of the Prevotellaceae family members in different analyses. Prevotella as a genus is one of the most abundant members in the gut microbiome of humans and several mammals. It is part of the abundant phylum Bacteroidetes, and its abundance is often inversely correlated with that of genus Bacteroides [127]; however, its reported biological effects are quite contra-dicting. Different members of family Prevotellaceae were identified throughout the analysis, but not all members or genera were found as statistically significantly different between the two time points or treatment groups. These variabilities can be resolved, in extended studies, by shotgun metagenomics or by long-read, full-length16S rRNA gene sequencing.
These types of limitations are common in multi-omics studies because data points are so numerous that statistical analyses often miss true-positive results. We consistently elected to err on the statistically conservative side to avoid false discovery, even at the expense of missing some true findings that remain to be resolved and confirmed in future studies.

Future Prospects
The findings of this study can be expanded in many directions. For example, testing whether pomegranate juice has a potential therapeutic effect can be achieved by repeating the experimental model with BaP administration for several weeks, followed by subsequent supplementation of pomegranate, after pathological damage has been confirmed (in a subset of animals) and BaP exposure has been discontinued.
In addition to testing the therapeutic potential of pomegranate juice, the same model can be expanded to test the intermediate and long-term effects of BaP, Pom, and PomBaP mixture in comparison to untreated animals, over several weeks or months.
Another way to expand the study is to use different fractions or components of pomegranate juice, or a combination thereof, to delineate the impact of each fraction and the potential synergistic activities of different bioactive components. Once pomegranate and its derivatives have been well tested, the model can test several other known or potential nutraceuticals.
Monitoring of the biotransformation of BaP and pomegranate in vivo via LC-MS or labeling is another direction for future studies.
From a translational perspective, the study is promising for the use of pomegranate juice, or its derivatives, as nutraceuticals. Two challenges need to be addressed in that case. First, the interactions of pomegranate with other drugs (reviewed in [128]) need to be taken into consideration. Second, the poor bioavailability of hydrolyzable tannins, which are key bioactive components of pomegranate juice, may be subject to formulation studies. Although the gut microbiota is pivotal in improving the bioavailability of hydrolyzable tannins, possible formulation alternatives include nanodelivery [129], micronization [130], and incorporation into bioactive film [131]. These alternatives may minimize microbiomedependent intra-and inter-individual variation in chemoprotection outcome.

Study Design
To assess the interaction between the pro-carcinogenic BaP, anti-carcinogenic pomegranate juice, and the rat gut microbiota, we randomly allocated 48 male Sprague-Dawley rats to four treatment groups: (i) an untreated group, only fed a standard AIN76 diet; (ii) a pomegranate group, fed the standard diet + 50 mg gallic acid equivalent (GAE)/kg/day pomegranate juice; (iii) a BaP group, fed the standard diet + 150 µg/kg/day BaP, dissolved in peanut oil [132]; and (iv) a pomegranate + BaP group, fed the standard diet + 50 mg GAE/kg/day pomegranate + 150 µg/kg/day BaP (dissolved in peanut oil). Thus, the 48 rats were assigned to four treatment groups, each represented by three cages and each housed by four rats.
All doses were calculated according to the average weights of rats within a given cage. Rats in all treatment groups received peanut oil either as a solvent for BaP, in groups iii and iv, or by oral gavage, in groups i and ii.

Animals
Seven-week-old male Sprague-Dawley rats (n = 48) were obtained from the modern veterinary office, Cairo, Egypt, and were housed in plastic-bottom cages (four per cage) with hardwood chip bedding and acclimatized for 10 days before the experiment started.
Rats were kept under standard conditions (temperature 25 ± 2 • C, relative humidity 60%) and a 12 h light/dark cycle. Food (AIN76 pellets) and water were given ad libitum.
The AIN76 diet was obtained from the Theodor Bilharz Research Institute, Cairo, Egypt. Pomegranate juice, peanut oil, and BaP were given to the rats daily by oral gavage for four weeks.
Animals were weighed weekly until the end of the experiment.

Preparation and Standardization of Pomegranate Juice
The entire edible part of the pomegranate fruit (Punica granatum L.), including the seeds (arils) and the attached sarcotestas, was thoroughly pressed in a Black & Decker Juice Extractor (Baltimore, MD, USA) to obtain concentrated pomegranate juice.
The juice was standardized by determination of total phenolics (TPC) according to the Folin-Ciocâlteu method [134]: A standard solution of gallic acid (10 mg/mL) was prepared by dissolving 1 g of gallic acid in 100 mL of methanol. A standard gallic acid curve was then constructed by preparing dilutions of 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.5 and 2 mg/mL in methanol from the standard solution (10 mg/mL) of gallic acid.
One hundred microliters of each dilution was mixed with 500 µL of distilled water and 100 µL of Folin-Ciocâlteu reagent. The mixture was then allowed to stand for 6 min. One milliliter of 7% sodium carbonate and 500 µL of distilled water were then added to the reaction mixture. After 90 min, the absorbance was spectrometrically measured at 760 nm [135]. The detailed characterization of individual phenolics in pomegranate juice by LC-MS has been previously described [136].

Sample Collection
At the start of the experiment, fecal samples were collected from each rat to evaluate basal glucuronidase enzyme activity. For sample collection, rats were placed in a clean, 70% alcohol-sterilized plastic box for no more than an hour. Stool was then aseptically collected by pre-sterilized forceps and placed in sterile 2 mL Eppendorf tubes. In some instances, sample collection required gently squeezing a rat's rectal region to excise the fecal pellet directly into sterile 2 mL Eppendorf tubes. Fecal samples were kept at −80 • C until use.
During the experiment, fecal samples were collected weekly and pooled from each cage into sterile 15 mL falcon tubes and kept at −80 • C until use.
At the end of the experiment, rats were sacrificed by cervical dislocation, and their colons were retrieved as per the guidelines of Ruehl-Fehlert et al. [137] and preserved in 10% formalin.

Beta-Glucuronidase Assay
To determine the β-glucuronidase enzymatic activity in rat feces, we prepared fecalase from the feces according to a modified version [122] of a previously published method [138]. Briefly, 50-70 mg fecal pellet was suspended in a potassium phosphate buffer (0.01 M, pH 7.4) and homogenized for 30 min. The fecal suspension was then centrifuged at 2000 rpm for 10 min, and the resulting supernatant was centrifuged again at 15,000 rpm for 20 min. The supernatant (fecalase) was then used for the enzyme activity assay.

Histopathological and Histochemical Examination
All histopathological/histochemical procedures and examinations were conducted at the Department of Cytology and Histology, Faculty of Veterinary Medicine, Cairo University (Giza, Egypt). The examiner was blinded for the sample groups.
Colonic tissue samples were flushed and fixed in 10% neutral buffered formalin for 72 h. Samples were trimmed and processed in serial grades of ethanol and cleared in xylene. The samples were then infiltrated and embedded into paraplast tissue-embedding media. Five-micrometer-thick tissue sections were cut by rotatory microtome for demonstration of intestinal wall layers in different samples.
Tissue sections were stained by hematoxylin and eosin (H&E, Tokyo, Japan) as a general morphological examination staining method and by Alcian blue (pH 2.5) for quantitative analysis of goblet cells' reactive acidic mucins, then they were examined microscopically in a blinded manner. All standard procedures for sample fixation and staining were conducted according to Culling [140].
For the histochemical analysis of reactive mucin, six random non-overlapping fields of colonic mucosa from each group were selected and scanned for the determination of area-based percentage of reactive mucins and goblet cells in Alcian-blue-stained (pH 2.5) tissue sections [141].
All light microscopic examinations and data were obtained by the Leica Application module for histological analysis attached to the full HD microscopic imaging system (Leica Microsystems GmbH, Wetzlar, Germany).

DNA Extraction and Quantification
Twenty-four fecal samples from different treatment groups (representing weeks 0 and 4) were subjected to DNA extraction. DNA was extracted by the QIAamp DNA Stool Minikit for microbial analysis (Qiagen, Hilden, Germany) and immediately stored at −20 • C.
The nucleic acid concentration and purity (260/280 and 260/230 ratios) of each sample were quantified in a nanodrop spectrophotometer. For DNA sequencing, the DNA was further quantified in a Qubit ® fluorometer (Life technologies, Carlsbad, CA, USA) by the Qubit dsDNA Hs assay kit (Life Technologies Corporation, Eugene, OR, USA).

16S rRNA Amplicon Sequencing
The concentrated DNA, of high quality and concentration, from each of the fecal samples, was sequenced at The Egyptian Center for Genome and Microbiome Research, Cairo, Egypt, in an iSeq™100 platform (Illumina, San Diego, CA, USA). The 2 × 150 bp paired-end protocol was followed, and an Illumina amplicon library was generated according to the Illumina16S rRNA library preparation protocol. The variable regions V3-V4 of bacterial 16S rRNA gene were amplified with degenerate PCR primers 515F

Microbiome Analysis
The raw sequence reads obtained from the iSeq instrument, after automated base calling by built-in BaseSpace software (Version 7.1.0) (Illumina, San Diego, CA, USA), were primarily analyzed by the Mothur software package v.1.36.1 [142] based on the MiSeq standard operating procedure (URL: https://mothur.org/wiki/miseq_sop, accessed on 25 April 2023) [143,144]. Short sequences, sequences that had unknown base pairs (sequences with N's), and those with at least one ambiguous base or more than eight homopolymers were removed.
Mothur taxonomic assignments were made against version 123 of the SILVA bacterial database [145] based on the Wang approach [146]. Sequences were then clustered into operational taxonomic units (OTUs) based on 97% sequence similarity cutoff, and chimeric sequences were identified and removed using UCHIME [147]. Of note, SILVA v. 123 uses traditional taxon names, some of which were updated in 2021 [148]. The use of traditional names allows comparison with current and past studies; however, we provide a list of renamed phyla (Table S2).
OTU abundance and consensus taxonomy output files were analyzed by Microbiome-Analyst (version 2.0)-a web-based tool for the statistical and visual analysis of microbiota data [65]. Reads were initially denoised by filters of a minimum read filtering of 20% prevalence with a count of 4. Sample reads were normalized and rarefied to the minimum library size, then data scaling (total sum scaling) was performed before the final differential abundance analysis. A low-variance filter was set at 5% for standard deviation. Alpha diversity profiling was conducted by the observed OTUs, Chao1 richness estimators, and Simpson and Shannon's indices, while statistical significance was tested by the Kruskal-Wallis method.
Beta diversity analysis was assessed by the permutational multivariate analysis of variance (PERMANOVA) statistical method, based on Bray-Curtis dissimilarity, and visualized by principal coordinate plots. Linear discriminant analysis effect size (LEfSe) was used to identify taxa with significant differences with default cutoffs. p values were adjusted by the false discovery rate (FDR) method [149].

Fecal Metabolite Profiling by Gas Chromatography Coupled with Mass Spectrometry Detection (GC-MS)
Seventy-two stool samples (triplicates of 24 samples) were processed for GC-MS metabolite analysis. For each sample, 750 µL of cold ethanol/20 mM phosphate extraction buffer (85:15 v/v) was added to 250 mg stool, and the mixture was vortexed for 2 min, homogenized for 2 min, and then centrifuged at 15,000 rpm at 4 • C for 15 min. The supernatant was then collected, and a 100 µL aliquot was mixed with 5 µL xylitol (1 mg/mL, as internal standard) and incubated with 200 µL of cold acetonitrile for 30 min at 4 • C, then centrifuged at 15,000 rpm for 10 min. The supernatant (100 µL) was collected and completely dried in a speed vacuum concentrator (Labconco Centrivap, Brermen, Germany).
Samples were then analyzed in a Shimadzu GC-17A gas chromatograph coupled with a Shimadzu QP5050A mass spectrometer. Silylated derivatives were separated on an Rtx-5MS (30 m length, 0.25 mm inner diameter, and 0.25 lm film) column. Samples were injected under 1:15 split-mode conditions: injector, 280 • C; column oven set at 80 • C for 2 min; and rate, 5 • C/min to 315 • C. Samples were kept at 315 • C for 12 min, with the helium carrier gas flow set at 1 mL/min. The transfer line and ion-source temperatures were set at 280 and 180 • C, respectively.

GC-MS Data Analysis
MS peak abundance of silylated metabolites were extracted with the MS Dial program, following the exact parameters settings for GC-MS [150]. The SIMCA-P version 14.1 software package (Umetrics, Umeå, Sweden) was used to export the aligned peak abundance data table, which was further exported for PCA and OPLS-DA [151]. All variables were mean-centered and scaled to Pareto variance.

Data Visualization and Statistical Analysis
Overall, normally distributed data with sample size greater than six were analyzed by parametric tests for significance (Student's t-test and ANOVA, for double or multiple variables, respectively), while data from small samples (six or under) or non-normally distributed data were analyzed by Wilcoxon and Kruskal-Wallis tests for double or multiple variables, respectively. For two-variable comparisons, paired tests were used for week 0 vs. week 4 comparisons, while unpaired tests were used to compare treatment groups. For multiple variable comparisons, post-hoc pairwise tests were performed with Bonferroni or FDR p value adjustment for smaller or larger sample sizes, respectively. Dunn test was used as an alternative to Bonferroni adjustment for non-parametric data. In some instances, both ANOVA and Kruskal-Wallis tests were applied when the latter was not applicable because of ties.
For microbiome analysis, the LEfSe tool was used to identify differences in relative abundance between different experimental groups, using an alpha value of 0.05, followed by the Kruskal-Wallis test with a threshold of 2.5 for logarithmic linear discriminant analysis (LDA) scores. PERMANOVA was used for beta diversity analysis, as detailed above.
The statistics were performed in version 4.1.2 of R [152] and visualized in RStudio [ [153], version 2022.02.0, Build 443] or MicrobiomeAnalyst (which uses R scripts). Correlations were identified by Spearman's rank correlation coefficient. A list of R packages used is provided in Table S3.

Sequence Deposition
All sequences were deposited to the NCBI sequence-reads archive (SRA), under project # PRJNA952221 and biosample numbers SAMN34069103-SAMN34069126.

Conclusions
In conclusion, this study allowed a head-to-head comparison of the effects of pomegranate and that of BaP on colonic tissues, mucin content, and on the general balance of the colonic microbiota and metabolome of adult Sprague-Dawley rats. The four weeks of the experiment, representing a substantial time in the life span of an adult rat, allowed the demonstration of animal age as a major predictor of microbiome and metabolome variations. In addition, the paired experimental design allowed controlling for inter-individual variations and cage effects by self-comparing each animal or cage in week 4 vs. week 0, thus leading to the resolution of treatment-derived variations. Finally, using various correlation and clustering analyses underscored microbiome-metabolome associations and highlighted microbiome co-abundance clusters and potential metabotypes that represent prospects for biomarker discovery. This work establishes a model for BaP-induced colon pathologies to assess the potential microbiota-mediated chemoprotective role of pomegranate juice (here) and other nutraceuticals and functional foods (in future studies).