Aryl Hydrocarbon Receptor (AhR) Activation by 2,3,7,8-Tetrachlorodibenzo-p-Dioxin (TCDD) Dose-Dependently Shifts the Gut Microbiome Consistent with the Progression of Non-Alcoholic Fatty Liver Disease

Gut dysbiosis with disrupted enterohepatic bile acid metabolism is commonly associated with non-alcoholic fatty liver disease (NAFLD) and recapitulated in a NAFLD-phenotype elicited by 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) in mice. TCDD induces hepatic fat accumulation and increases levels of secondary bile acids, including taurolithocholic acid and deoxycholic acid (microbial modified bile acids involved in host bile acid regulation signaling pathways). To investigate the effects of TCDD on the gut microbiota, the cecum contents of male C57BL/6 mice orally gavaged with sesame oil vehicle or 0.3, 3, or 30 µg/kg TCDD were examined using shotgun metagenomic sequencing. Taxonomic analysis identified dose-dependent increases in Lactobacillus species (i.e., Lactobacillus reuteri). Increased species were also associated with dose-dependent increases in bile salt hydrolase sequences, responsible for deconjugation reactions in secondary bile acid metabolism. Increased L. reuteri levels were further associated with mevalonate-dependent isopentenyl diphosphate (IPP) biosynthesis and o-succinylbenzoate synthase, a menaquinone biosynthesis associated gene. Analysis of the gut microbiomes from cirrhosis patients identified an increased abundance of genes from the mevalonate-dependent IPP biosynthesis as well as several other menaquinone biosynthesis genes, including o-succinylbenzoate synthase. These results extend the association of lactobacilli with the AhR/intestinal axis in NAFLD progression and highlight the similarities between TCDD-elicited phenotypes in mice to human NAFLD.


Introduction
Non-alcoholic fatty liver disease (NAFLD) is estimated to affect~25% of the global population and is defined as a spectrum of progressive pathologies, including steatosis, immune cell infiltration/inflammation, fibrosis, and cirrhosis. It is associated with increased risk for hepatocellular carcinoma and is the second leading cause of liver transplants in the USA [1]. Other pathologies, including obesity, type 2 diabetes (T2D), and coronary heart disease, demonstrate a high co-occurrence with NAFLD, e.g.,~40-70% in T2D patients and 90% in obese patients [2]. A multi-hit hypothesis for NAFLD proposes several contributing factors to development and progression, including disruptions in the immune system, adipose tissue metabolism, and the gut microbiome [3].
Previous work demonstrated that the serum levels of LCA and DCA increased following TCDD treatment suggesting enrichment for the microbial bile acid metabolism [9]. To further explore dose-dependent disruptions in the gut microbiome and microbial metabolism relevant to the progression of NAFLD-like pathologies, shotgun metagenomic analysis was used to examine the dose dependent taxonomic and metabolic disruptions elicited by TCDD.

TCDD-Elicited Toxicity Enriched for Lactobacillus Species
Taxonomic analysis identified significant dose-dependent population shifts among caecum microbiota in response to TCDD. While no significance was observed between treatment groups at the phylum level, a decreasing trend was observed for Bacteroidetes concurrent with increasing trends in Firmicutes abundance ( Figure 1A). In regard to the host, some secondary bile acids, e.g., lithocholic acid [LCA] and deoxycholic acid [DCA]), exhibit a high affinity for the farnesoid x receptor (FXR) and G proteincoupled bile acid receptor (TGR5, a.k.a., GPBAR1), which regulate glucose, lipid, and bile acid homeostasis [34][35][36]. In human NAFLD, the secondary bile acid metabolism is disrupted with bile acid analogs that target the FXR and TGR5 signaling pathways under development for the treatment of liver disease [22,31,37]. Previous work demonstrated that the serum levels of LCA and DCA increased following TCDD treatment suggesting enrichment for the microbial bile acid metabolism [9]. To further explore dose-dependent disruptions in the gut microbiome and microbial metabolism relevant to the progression of NAFLD-like pathologies, shotgun metagenomic analysis was used to examine the dose dependent taxonomic and metabolic disruptions elicited by TCDD.

TCDD-Elicited Toxicity Enriched for Lactobacillus Species
Taxonomic analysis identified significant dose-dependent population shifts among caecum microbiota in response to TCDD. While no significance was observed between treatment groups at the phylum level, a decreasing trend was observed for Bacteroidetes concurrent with increasing trends in Firmicutes abundance ( Figure 1A). Taxa abundance were assessed in metagenomic cecum samples from male C57BL/6 mice following oral gavage with sesame oil vehicle or 0.3, 3, or 30 µg/kg TCDD every 4 days for 28 days (n = 3). Significant shifts in relative abundances of taxa are presented at the (A) phylum, (B) genus, (C) and species levels. Significance is denoted with an asterisk (*; adjusted p-value < 0.1). Figure 1. TCDD enriched Lactobacillus species in the cecum microbiota. Taxa abundance were assessed in metagenomic cecum samples from male C57BL/6 mice following oral gavage with sesame oil vehicle or 0.3, 3, or 30 µg/kg TCDD every 4 days for 28 days (n = 3). Significant shifts in relative abundances of taxa are presented at the (A) phylum, (B) genus, (C) and species levels. Significance is denoted with an asterisk (*; adjusted p-value < 0.1).
At the genus level, Turicibacter was enriched by TCDD while the genus Lactobacillus trended towards enrichment ( Figure 1B). Interestingly, at the species level, 10 out of 13 enriched species were from the Lactobacillus genus (e.g., L. reuteri and Lactobacillus sp. ASF360) as well as Turicibacter sanguinis. Conversely, the most abundant Lactobacillus species in vehicle treated mice, Lactobacillus murinus, trended towards a dose-dependent decrease ( Figure 1C). The major changes in taxa were observed in the 30 µg/kg TCDD treatment group.

Bile Salt Hydrolase (Bsh) Levels Correlated with Significantly Enriched Species
Many Lactobacillus species deconjugate primary conjugated bile acids mediated by bile salt hydrolases (BSH), imparting bile acid tolerance [38]. To further investigate the effect of TCDD on bile acid metabolism, bsh sequences were annotated and quantified within metagenomic samples. Annotations to bsh were increased by TCDD and associated with enriched species, including L. reuteri and T. sanguinis ( Figures 1C and 2A, and Table S1). At the genus level, Turicibacter was enriched by TCDD while the genus Lactobacillus trended towards enrichment ( Figure 1B). Interestingly, at the species level, 10 out of 13 enriched species were from the Lactobacillus genus (e.g., L. reuteri and Lactobacillus sp. ASF360) as well as Turicibacter sanguinis. Conversely, the most abundant Lactobacillus species in vehicle treated mice, Lactobacillus murinus, trended towards a dose-dependent decrease ( Figure 1C). The major changes in taxa were observed in the 30 µg/kg TCDD treatment group.

Bile Salt Hydrolase (Bsh) Levels Correlated with Significantly Enriched Species
Many Lactobacillus species deconjugate primary conjugated bile acids mediated by bile salt hydrolases (BSH), imparting bile acid tolerance [38]. To further investigate the effect of TCDD on bile acid metabolism, bsh sequences were annotated and quantified within metagenomic samples. Annotations to bsh were increased by TCDD and associated with enriched species, including L. reuteri and T. sanguinis ( Figures 1C and 2A, and Table  S1). . The presence of bsh gene sequences were assessed in metagenomic caecum samples from male C57BL/6 mice following oral gavage with sesame oil vehicle or 0.3, 3, or 30 µg/kg TCDD every 4 days for 28 days using three independent cohorts (n = 3). (A) The presence (green boxes) or absence of bsh sequences detected in any of the metagenomic samples (n = 3) are denoted within the respective treatment groups. Significant increases (*) or decreases (@) in normalized bsh abundances (adj. p < 0.1) are denoted. Also denoted is significantly increased species (#) determined by taxonomic analysis that . The presence of bsh gene sequences were assessed in metagenomic caecum samples from male C57BL/6 mice following oral gavage with sesame oil vehicle or 0.3, 3, or 30 µg/kg TCDD every 4 days for 28 days using three independent cohorts (n = 3). (A) The presence (green boxes) or absence of bsh sequences detected in any of the metagenomic samples (n = 3) are denoted within the respective treatment groups. Significant increases (*) or decreases (@) in normalized bsh abundances (adj. p < 0.1) are denoted. Also denoted is significantly increased species (#) determined by taxonomic analysis that corresponded with respective RefSeq species bsh annotations. Significance was determined by Maaslin2 R package. (B) Volcano plot displaying log2(fold-changes) in relative abundance of species between vehicle and 30 µg/kg TCDD treatment groups versus -log(adjusted p-values [adj. p]). Red dots denotes bsh sequences detected in 30 µg/kg TCDD treatment group. Significance was determined by the DeSeq2 R package comparing only vehicle and 30 µg/kg TCDD groups. Red dashed lines are reference to −log(0.05) value for the y-axis and −1 and 1 for the x-axis.
Conversely, L. murinus associated bsh annotations exhibited a dose-dependent decrease consistent with decreasing trends in taxonomic abundance. Although not reaching significance, many bsh sequences were also associated with unclassified Lachnospiraceae species, including Lachnospiraceae bacterium A4, a community member reaching 5-23% relative abundance in the cecum metagenomic samples (Figure 2A). In contrast, Lactobacillus gasseri was enriched but no bsh sequences were identified ( Figure 2B). To summarize, the top enriched species were also associated with increased abundances in bsh levels in the cecum.

TCDD Enriched for Mevalonate-Dependent Isoprenoid Biosynthesis
To investigate other metabolic pathways imparting competitive advantages to TCDDelicited gut environmental stresses, functional gene annotations associated with L. reuteri, the highest enriched species, were assessed. Among enriched uniref90 annotations in the cecum metagenomic dataset was the aromatic amino acid aminotransferase (UniRef90_ A0A2S1ENB9) also classified to L. reuteri (Table S2). Aromatic amino acid aminotransferase produces a tryptophan metabolite, indole-3-aldehyde, a known AhR ligand reported to induce IL-22 in vivo [26]. Among 39 enzyme commission (EC) annotations that were enriched and associated with L. reuteri were several annotated to the isoprenoid biosynthesis pathway (Figure 3, Figures S1 and S3).
Bacteria biosynthesize the isoprenoid, isopentenyl diphosphate (IPP), either through the mevalonate-dependent pathway, which is also found in mammals, or the methylerythritol phosphate (MEP)-pathway. Both L. reuteri and Lactobacillus johnsonii were the major contributors to mevalonate-dependent IPP biosynthesis pathway enrichment with almost all genes in the pathway increased by TCDD; four out of six of the genes significantly increased by TCDD (Figure 3 and Figure S1). Gene enrichment in the alternative MEPpathway were unchanged by TCDD. For L. murinus, only two EC annotations (EC 2.7.1.148, 4-diphosphocytidyl-2-C-methyl-D-erythritol (CDP-ME) kinase, and EC 5.3.3.2, isopentenyldiphosphate Delta-isomerase) were identified in the MEP pathway also found in L. reuteri ( Figure S2). HUMAnN 3.0 analysis of a published metagenomics dataset of fecal samples from human cirrhotic patients (https://www.ncbi.nlm.nih.gov/bioproject/PRJEB6337/, accessed on 25 March 2021) [39] revealed strikingly similar results to our caecum samples from TCDD treated mice. Specifically, increased gene abundance associated with the mevalonatedependent pathways was also evident in patients with compensated and decompensated liver disease (Figure 4). in the cecum metagenomic dataset was the aromatic amino acid aminotransferase (UniRef90_A0A2S1ENB9) also classified to L. reuteri (Table S2). Aromatic amino acid aminotransferase produces a tryptophan metabolite, indole-3-aldehyde, a known AhR ligand reported to induce IL-22 in vivo [26]. Among 39 enzyme commission (EC) annotations that were enriched and associated with L. reuteri were several annotated to the isoprenoid biosynthesis pathway (Figures 3, S1, and S3).  Compensated cirrhosis is defined as no decrease in liver function while decompensated cirrhosis exhibits decreased liver function. Among decompensated patients with cirrhosis, the mevalonate dependent IPP pathway was increased in 7 out of 8 EC numbers required for de novo IPP biosynthesis ( Figure 4). Taxa annotated to genes in the pathway exhibited a wide variety in genera for each EC number in human samples compared to murine cecum samples from this study ( Figure S3). Taxonomy classified to a majority of the mevalonate-dependent genes were from the Lactobacillaceae family, including Enterococcus, Lactobacillus, Streptococcus genera ( Figure S3 and Table S5). Lactobacillus and Streptococcus species, including L. reuteri and Streptococcus anginosus, a known pathogen in liver abscesses [40], were among species classified to the pathway (Tables S4 and S5). HUMAnN 3.0 analysis of a published metagenomics dataset of fecal samples from human cirrhotic patients (https://www.ncbi.nlm.nih.gov/bioproject/PRJEB6337/, accessed on 25 March 2021) [39] revealed strikingly similar results to our caecum samples from TCDD treated mice. Specifically, increased gene abundance associated with the mevalonate-dependent pathways was also evident in patients with compensated and decompensated liver disease ( Figure 4).

Phenotypes and Gut Microbiomes of Cirrhosis Patients
In polyprenol diphosphate biosynthesis, IPP is recursively added to geranyl diphosphate (GPP) or farnesyl diphosphate (FPP) for polyprenol biosynthesis used in vitamin K2 (a.k.a., menaquinone) and peptidoglycan biosynthesis [41,42]. TCDD enriched for heptaprenyl diphosphate synthase (EC 2.5.1.30) with major contributions from L. reuteri and L. johnsonii ( Figure 5).  As bacterial cell wall restructuring has been reported in response to bile acids and different levels of isoprenoid biosynthesis pathways were identified, peptidoglycan biosynthesis was also assessed [43]. Most genes encoding enzymes required for peptidoglycan biosynthesis were present in the metagenomic dataset ( Figure 6A) with no changes observed following TCDD treatment.
As bacterial cell wall restructuring has been reported in response to bile acids and different levels of isoprenoid biosynthesis pathways were identified, peptidoglycan biosynthesis was also assessed [43]. Most genes encoding enzymes required for peptidoglycan biosynthesis were present in the metagenomic dataset ( Figure 6A) with no changes observed following TCDD treatment.   In the cirrhosis samples, several EC numbers representing the initial menaquinon biosynthesis steps were also increased in compensated and decompensated patients (Fig  ure 8, steps 1  In the mouse cecum dataset, species contributing to o-succinylbenzoate menaquinone biosynthesis pathway included Escherichia coli, several Bacteroides (e.g., Bacteroides vulgatus and Bacteroides caecimuris), and Lactobacillus species (e.g., L. reuteri) ( Figure 7B and Table S6). No one species was annotated to the entire set of enzymes needed for de novo biosynthesis from chorismate, however B. vulgatus was annotated for 6 out of 9 genes in the pathway (Table S6). O-Succinylbenzoate synthase ( Figure 6A, EC 4.2.1.113, step 4) was increased by 30 µg/kg TCDD with L. reuteri being the major contributor to relative abundance (Figure 7, step 4).
Lactobacillus species annotated to menaquinone biosynthesis included L. reuteri, L. murinus, and L. johnsonii. Among annotated menaquinone biosynthesis EC numbers, L. reuteri was among the identified Lactobacillus species that had the highest relative abundance and most menaquinone EC annotations ( Figure 7B). L. reuteri also had annotations in samples for EC numbers involved in the final steps of the shikimate pathway responsible for chorismate biosynthesis (Table S7).
In the cirrhosis samples, several EC numbers representing the initial menaquinone biosynthesis steps were also increased in compensated and decompensated patients (Figure 8, steps 1,3-5), including o-succinylbenzoate synthase (Figure 8, EC 4.2.1.113, step 4). However, L. reuteri was not among species classified to this EC number. Species classified to all EC numbers comprising the complete pathway included E. coli, and Klebsiella species, such as K. pneumoniae and Citrobacter species. L. reuteri, were not annotated to any menaquinone biosynthesis genes in healthy or compensated patients, but several EC numbers in the decompensated group (EC 6.2.1.26, 4.1.3.6, and 2.1.1.163), which are involved in later stages of menaquinone biosynthesis (Table S8).

Discussion
Previous studies have reported that TCDD elicited NAFLD-like pathologies, dysregulated bile acid metabolism and gut microbiome dysbiosis [9,11,12,28,30]. This study fur- However, L. reuteri was not among species classified to this EC number. Species classified to all EC numbers comprising the complete pathway included E. coli, and Klebsiella species, such as K. pneumoniae and Citrobacter species. L. reuteri, were not annotated to any menaquinone biosynthesis genes in healthy or compensated patients, but several EC numbers in the decompensated group (EC 6.2.1.26, 4.1.3.6, and 2.1.1.163), which are involved in later stages of menaquinone biosynthesis (Table S8).

Discussion
Previous studies have reported that TCDD elicited NAFLD-like pathologies, dysregulated bile acid metabolism and gut microbiome dysbiosis [9,11,12,28,30]. This study further elucidated the shifts in the gut microbiota associated with TCDD treatment using shotgun metagenomic sequencing. We show that TCDD dose-dependently shifted the gut microbiota composition by enriching for Lactobacillus species, consistent with hepatic disruption of host and microbial bile acid metabolism. In addition, TCDD enriched for genes involved in mevalonate dependent isoprenoid precursor biosynthesis and menaquinone biosynthesis, crucial for microbial cell growth and survival. Over-representation of these microbial associated pathways were also identified in human cirrhosis stool metagenomics datasets.
We observed a species-specific increase of L. reuteri with a concurrent decrease in L. murinus suggesting shifts in Lactobacillus composition at the species and/or strain levels. Further, decreased abundance of L. murinus has been reported in human NAFLD [53]. Other taxa enriched following treatment included Turicibacter sanguinis, an anaerobic grampositive bacillus commonly found in animals, including humans [54]. Interestingly, T. sanguinis has been shown to deconjugate bile acids and metabolize serotonin affecting lipid and steroid metabolism [54,55].
Quantitative trait locus analysis correlated T. sanguinis abundance with cholic acid levels and expression of the intestinal bile acid transporter Slc10a2 [54]. Both cholic acid levels and Slc10a2 expression are dose-dependently increased by TCDD [9]. Consequently, the dose-dependent taxonomic shift in Lactobacillus and Turicibacter species known to deconjugate conjugated bile acids is consistent with increased levels of secondary bile acids following TCDD treatment.
Some host relevant intestinal health and homeostatic effects can be attributed to Lactobacillus species mediated by bile salt hydrolases (BSHs), which are responsible for deconjugation reactions, the gateway step for conversion of conjugated primary bile acid to secondary bile acids [56]. A majority of Lactobacillus species possess BSHs, often containing multiple different gene copies within their genome, some with different bile acid substrate preferences [33,38].
However, the presence of bsh sequences does not simply infer bile acid tolerance as growth inhibition and reduced fitness is also possible depending on the conjugated or deconjugated bile acids present and/or BSH specificity [33,38,57]. For example, L. gasseri bsh knockout mutants exhibit increased fitness compared to wild type strains [38]. Interestingly, L. gasseri bsh sequences were not identified despite increased L. gasseri abundance following TCDD treatment. Our bsh analysis also found TCDD enriched Lactobacillus-associated sequences that may impart bile acid tolerance.
For example, the bsh sequence enriched by TCDD annotated to L. johnsonii (RefSeq ID: EGP12391) (Table S3) exhibited higher substrate specificity for glycine over taurine conjugated bile acids [58]. In a companion study using the dose response and treatment regimen, Fader et al. reported TCDD increased serum DCA levels~80 fold, with only ã two-fold increase in serum GDCA levels [9]. In contrast, hepatic taurolithocholic acid (TLCA) levels were increased~233 fold while serum lithocholic acid increased only four fold following TCDD treatment. Moreover, glycine conjugated bile acids, including GDCA, are more toxic towards Lactobacillus species than taurine conjugated bile acids [33,59,60]. Increased levels of BSH with a substrate preference for glycine conjugated bile acid may partially explain select Lactobacillus species enrichment. Further, both TLCA and DCA are potent FXR and GPBAR1 agonists that regulate the lipid, glucose, and bile acid metabolism [61,62]. Consequently, shifts in microbial secondary bile acids by Lactobacillus species may play a role in TCDD elicited gut dysbiosis impacting host regulation of energy homeostasis.
Coincident with increased levels of bsh was the dose-dependent increase in genes from the mevalonate-dependent isoprenoid biosynthesis, the pathway also used in mammals for cholesterol biosynthesis. The MEP pathway is the predominant isoprenoid biosynthesis pathway among gut microbiota while the mevalonate-dependent pathway is only found in select bacteria, including Lactobacillus and Streptococcus species [63]. The output from either pathway is farnesyl diphosphate (FPP) and geranyl diphosphate (GPP), substrates required for polyprenol biosynthesis used in menaquinone and cell wall biosynthesis.
Menaquinones are utilized by bacteria for anaerobic/aerobic respiration, providing antioxidant activity with menaquinone supplementation affecting the gut microbiome [64]. In the context of L. reuteri, we observed genes annotated to the shikimate pathway, which is responsible for chorismate biosynthesis, a precursor for aromatic amino acids and the naphthoquinone head group of menaquinone, as well as genes involved in de novo menaquinone biosynthesis. While the complete biosynthesis pathway was not present in L. reuteri, it is consistent with other metagenomic reports of incomplete menaquinone biosynthesis pathways in gut Lactobacillus species [44].
It has been proposed that Lactobacillus species may participate in later menaquinone biosynthesis steps through the uptake of intermediates, such as o-succinylbenzoate from other bacteria or dietary sources [44]. In addition, the ability to utilize menaquinones for respiration is typically not associated with Lactobacillus species. However, some lactic acid bacteria, including L. reuteri strains, demonstrate the ability to respire when menaquinone and heme are supplemented [65,66].
Metagenomic analysis also identified the mevalonate-dependent pathway enrichment in fecal samples from patients with cirrhosis. The mevalonate-dependent pathway is reported to be increased in fibrosis patients with autoimmune pathologies [67]. Isoprenoid biosynthesis pathways are also elevated in the lung microbiome of cystic fibrosis patients, with the MEP pathway enriched rather than the mevalonate route [68]. The association between fibrosis and isoprenoid biosynthesis enrichment warrants further investigation in the context of potential mechanisms contributing to bacterial fitness and/or fibrosis.
Increased abundance of the mevalonate-dependent biosynthesis pathway could also be a biomarker of Lactobacillus and Streptococcus proliferation that is often associated with non-alcoholic steatohepatitis (NASH)/fibrosis [21,45]. We identified enrichment of the mevalonate-dependent pathway in both mouse and human microbiomes, whereas the complete pathway was primarily annotated to Streptococcus and Lactobacillus species (Table S7).
Other factors, such as simvastatin and proton pump inhibitors (PPI), that are commonly prescribed for NAFLD patients may also impact these microbial pathways. Simvastatin, which is primarily excreted in the feces [69], has been reported to reduce bacterial growth by directly inhibiting bacterial HMG-CoA synthesis while PPIs inhibit Streptococcus species growth [70][71][72]. These microbiome-drug interactions highlight off target effects that should be considered when investigating novel NAFLD treatments, such as new drug development and/or probiotic interventions.
In addition to increased mevalonate-dependent isoprenoid biosynthesis genes in cirrhotic patients, menaquinone biosynthesis gene abundance was also increased. This suggests taxa with the ability to produce menaquinone may have a competitive advantage when intestinal environmental conditions shift during disease progression. In cirrhosis patients, E. coli and B. vulgatus were associated with genes providing a majority of the menaquinone biosynthesis capacity. These species are also increased in human NAFLD [73].
Similar to the results in mice exposed to TCDD, L. reuteri was associated with several menaquinone biosynthesis genes and only detected in decompensated cirrhosis patients but lacked the complete pathway (Table S10). In cirrhosis patients, it is unclear whether L. reuteri is participating in menaquinone metabolism and/or benefiting from increased abundance of species, like E. coli that are capable of producing menaquinones.
This study was designed to account for factors affecting gut microbiota analysis, including cage effects, coprophagia, and circadian rhythms [74]. Significant shifts in taxa were observed in Lactobacillus species. However, the small group size (n = 3) following adjustment for multiple testing lacked sufficient power to confirm more subtle shifts, such as the two fold enrichment of Lachnospiraceae A4, an abundant community member associated with bsh sequences.
Samples were also collected within the same Zeitgeber period to account for possible variations in relative microbiota levels due to circadian rhythm/diurnal regulation. In fact, L. reuteri is one gut microbiome member demonstrating changes in relative abundance in human samples due to circadian/diurnal regulation [75]. TCDD disrupted diurnal regulation of hepatic gene expression, including bile acid biosynthesis genes, which may contribute to L. reuteri enrichment [76].
Using this study design, TCDD also disrupted the bile acid metabolism and enterohepatic circulation with increased hepatic and serum total bile acids and secondary bile acid DCA [9]. Likewise, increased serum bile acid levels, including DCA, have been reported in patients with steatohepatitis and fibrosis [79][80][81][82]. Furthermore, NAFLD patients with increased bile acid levels have increased levels of bacterial genes from the bai operon associated with 7α-dehydroxylation of bile acids leading to the production of DCA from cholic acid [83].
Although the consequences of TCDD-elicited immune system effects on the gut microbiome were not assessed in this study, it is most likely a factor impacting L. reuteri enrichment. TCDD causes macrophage and dendritic cell migration out of the lamina propria with increased accumulation in the liver, possibly exacerbating hepatic inflammation and affecting intestinal immune responses [14]. The ability of L. reuteri to produce AhR ligands, upregulate IL-22, and associate with the mucosa and Peyer's patches provides geographical proximity for immune/microbiome crosstalk mediated by the AhR [26,84,85].
In addition to immune cell regulation, TCDD increased bone formation and decreased bone marrow adiposity [86]. Interestingly, L. reuteri supplementation also increased bone density, but only when mice were induced towards an inflammatory state [87]. Overall, the dose-dependent increase in L. reuteri levels is consistent with increased bile acid levels, disruption of circadian/diurnal regulation and increased bone density [9,75,86,87].
In summary, Lactobacillus species were dose-dependently increased following AhR activation by TCDD concurrent with the increase in bsh genes and increased primary and secondary bile acids. Specifically, L. reuteri, a keystone gut microbiome species is involved in the microbial metabolism of bile acids and AhR ligands. The large and uniform enrichment of L. reuteri in this study also suggests environmental pressures, such as increased levels of bile acids and antimicrobial peptides elicited by AhR activation, may provide a complementary niche for L. reuteri that possess a gene repertoire not found in the closely related L. murinus.
We also provide evidence regarding how the L. reuteri metabolism could impact the AhR, FXR, and GPBAR1 signaling pathways, placing L. reuteri at the crossroads of bacterial/host interactions affecting glucose, bile acid, and immune regulation. Whether these microbial shifts in metabolism are adaptive and limit the intensity of adverse consequences or exacerbate steatosis to steatohepatitis with fibrosis progression warrants further investigation.
The rodent diet is a fixed formula complete diet with an energy density of 3.0 kcal/g and a nutrient ingredient composition, including 22% protein, 5.5% fat, and 40.6% carbohydrate. Mice (PND29) were orally gavaged at the beginning of the light cycle (between Zeitgeber time 0-3) with 0.1 mL sesame oil vehicle (Sigma-Aldrich, St. Louis, MO, USA) or 0.3, 3 and 30 µg/kg body weight TCDD (AccuStandard, New Haven, CT, USA) every 4 days for 7 total exposures (n = 3 per treatment group). The study was conducted in three cohorts with mice housed separately among treatment groups for a total of 9 mice per treatment group.
In each cohort, three mice were housed per treatment group, and one mouse was randomly selected from each treatment group per cohort (n = 3 per treatment group for the metagenomic analysis) to account for coprophagia and ensure reproducibility. The first gavage was administered on day 0 of the study. On day 28, vehicle-and TCDD-treated mice (fasted for 6 h with access to water) were weighed and euthanized by CO 2 inhalation between Zeitgeber time 0-3. Upon collection, cecums were immediately flash frozen in liquid nitrogen and stored at −80 • C until analysis. All animal handling procedures were performed with the approval of the Michigan State University (MSU) Institutional Animal Care and Use Committee.

Metagenomic Taxonomic Analysis
Kaiju was used for taxonomic analysis of mouse cecum metagenomic dataset. The reference database used was the progenomes database downloaded from the kaiju webserver (https://kaiju.binf.ku.dk/database/kaiju_db_progenomes_2020-05-25.tgz, accessed on 15 October 2020). Multivariate association between dose and taxonomy relative abundances used Maaslin2 (https://github.com/biobakery/Maaslin2, accessed on 15 October 2020) [91] with the following default settings used: normalization (total sum scaling), analysis method (general linear model), and Benjamini-Hochberg multiple test correction. Adjusted p-values for Maaslin2 analysis used dose (sesame oil vehicle (0), 0.3, 3, or 30 µg/kg TCDD) as the fixed effect, which was treated as continuous variable and the vehicle set for reference. For comparison of taxonomy between vehicle and 30 µg/kg TCDD treatment groups, DeSeq2 was used to determine adjusted p-values using default settings [92].

Metagenomic Functional Analysis
The HUMAnN 3.0 bioinformatic pipeline [93] was used with default settings to classify reads to UniRef90 protein identifications using UniProt's UniRef90 protein data base (January 2019, accessed on 15 January 2021). Reads aligned to UniRef90 identifications were mapped to enzyme commission (EC) number entries using the human_regroup_table tool. Read abundance was normalized to gene copies per million reads (CPM) using the human_renorm_table tool. Multivariate association between dose and enzyme commission number relative abundance used Maaslin2 with same settings used for taxonomy analysis.
Xander (a gene-targeted assembler, https://github.com/rdpstaff/Xander_assembler, accessed on 5 February 2021) was used to annotate and quantify bile salt hydrolase sequences with the following settings: k-mer size = 45, filter size = 35, minimum assembled contig bit score = 50, and minimum assembled protein contigs = 100 [94]. Reference DNA and protein bsh sequences used for Xander were downloaded from FunGenes Gene Repository and are listed in supplementary material (Tables S9 and S10) [95]. For RefSeq bsh sequence analysis, relative abundance was determined by normalizing to total abundance of rplB sequences also determined by Xander per sample. Significance was determined with Maaslin2 with the same settings used for taxonomy analysis.
Cirrhosis was diagnosed by either biopsy, clinical evidence of decompensation, or other metrics, including radiological evidence of liver nodularity and intra-abdominal varices in a patient with chronic liver disease [39]. The subclassification was used as fixed effect for analysis with healthy as the reference category. Again, Maaslin2 was used with settings used for mouse functional analysis with diagnosis as a fixed effect with healthy diagnosis as reference to determine adjusted p-values for compensated and decompensated patient designations.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the human metagenomics study [39].
Data Availability Statement: Quality filtered mouse metagenomic data from this study can be found at the NCBI Sequence Read Archive (SRA) (https://www.ncbi.nlm.nih.gov/bioproject/ PR-JNA719224/, accessed on 10 August 2021) under the accession ID PRJNA719224. The human fecal metagenomic data presented in this study is openly available and can be found at the NCBI SRA under the accession number PRJEB6337 (https://www.ncbi.nlm.nih.gov/bioproject/PRJEB6337/, accessed on 25 March 2021) [39].