Effect of berberine on global bile acid metabolome. Since previous reports have indicated that BBR affects hepatic lipids and cholesterol by increasing bile acid excretion into the large intestine , the global bile acid metabolome was examined in control and BBR treated gnotobiotic mice. Total liver bile acid concentrations were not significantly different between control and BBR treated mice (P = 0.5283) (Fig. 1A). Significant compositional differences between bile acids in BBR treated and control liver and serum were not observed; however, microbial secondary bile acid products such as deoxycholic acid (DCA), taurodeoxycholic acid (TDCA), and taurolithocholic acid (TLCA) were observed, indicating functional bile acid metabolism by the B4PC2 consortium in the GI tract (Figure S1-S3). By contrast, a significant increase in total cecal bile acids (4.57 ± 1.42 µmol/g vs. 1.29 ± 0.106 µmol/g; P < 0.05) (Fig. 1B), and cecal bile acid composition was observed after BBR treatment relative to the control (Fig. 1C & S4). These disparate responses to BBR treatment observed in liver, serum, and cecum suggest that BBR increases fecal loss of bile acids with concomitant increased synthesis of bile acids in order to maintain baseline liver bile acid concentrations.
Functional bile acid metabolism in the cecum by the human gut B4PC2 community was evident in both control and BBR treatment groups. Deconjugation of taurine-conjugated bile acids was nearly complete. DCA and LCA, the end-products of the bile acid 7α-dehydroxylation pathway encoded by C. hylemonae and C. hiranonis , were major metabolites in the cecum (Figure S4). Conversion of CDCA to α- and β-muricholic acid (MCA) was observed (Figure S4); however, little conversion of MCA to murideoxycholic acid (MDCA) and ω-MCA was detected. This confirmed the mice in this study housed a microbial consortium of human bacteria, as human mixed fecal bacteria and bile acid 7α-dehydroxylating clostridia appear to be unable to metabolize MCA to MDCA and ω-MCA [11, 13–15]. Additionally, MDCA has been shown to be generated by the host .
Cecal composition of the B4PC2 human microbial consortium in control and berberine treated gnotobiotic mice. To determine whether the B4PC2 consortium was established in both control and BBR treated mice a microbial analysis of cecal content was performed. Sequencing of 16S rDNA resulted in 53,456 ± 3,743 reads for control ceca and 53,529 ± 3,441 reads for BBR treated ceca. Overall diversity of both control and BBR mice were consistent with the inoculated consortium indicating successful maintenance of the germ-free environment (Figure S5). To examine whether the composition of this bile-tolerant community is altered by BBR treatment, diversity analyses were performed on 23,900 rarefied reads from the 16S rDNA dataset. Non-metric multidimensional scaling (NMDS) and Analysis of Similarities (ANOSIM) tests of differences in beta-diversity resulted in R = 0.141 and P = 0.123 (999 permutations) indicating diversity between samples was not significant (Figure S5). Shannon index (alpha diversity) was not significantly different between groups (P = 0.30; Mann-Whitney test) (Figure S6). We next performed network correlation analysis on bile acids in the cecum, serum, and liver with abundances of B4PC2 consortium members. Substantial network topological changes between bile acids, and microbial taxa were observed between control and BBR networks (Fig. 2; Supplementary Dataset 1). These analyses indicate that the B4PC2 consortium was similarly established in control and BBR treated gnotobiotic mice, suggesting that BBR mechanisms of action are not related to composition changes in this bile acid metabolizing microbial community.
Direct effects of berberine and bile acid concentrations on gene expression by the B4PC2 consortium. Given observed topological changes between the bile acid metabolome and the B4PC2 consortium, RNA-Seq was performed on cecal content collected from control and BBR treated mice. Network correlation analyses between transcriptomic and metabolomic data were created for each member of the B4PC2 consortium in order to decern whether differences in gene expression are in response to bile acid concentration or direct effect of BBR treatment. Results for each B4PC2 member are highlighted in the following sections.
Bilophila wadsworthia. Berberine treatment resulted in significant differential expression of 123 genes (74 up-regulated; 49 down-regulated) (Fig. 3A; Supplementary Dataset 2–4). Transcriptome data indicate that in presence of BBR, B. wadsworthia imports bacterial membrane lipids (phosphatidylethanolamine) (LadL; 3.29 log2FC, P = 0.01), degrades ethanolamine to acetaldehyde and ammonia (ethanolamine ammonia-lyase), and converts acetaldehyde to ethanol (adh; 3.29 log2FC, P = 2.46E-3). High relative expression of group 1b NiFeSe hydrogenase was observed (3.56 log2FC, P = 3.66E-4), which is involved in liberation of electrons for formate, sulfite, and nitrate respiration. Pyruvate-formate lyase and formate dehydrogenase were highly expressed in Bilophila during BBR treatment. However, the gene most highly-expressed was nitrate reductase γ-subunit (8.04 log2FC, P = 7.07E-06). The other two most highly expressed genes were citric acid cycle enzyme malate dehydrogenase (5.43 log2FC, P = 1.04E-5) and citrate transporter (5.11 log2FC, P = 1.04E-4). Additional citric acid cycle genes and respiratory complex genes are significantly up-regulated by BBR (Fig. 3A; Supplementary Dataset 2–4). In addition, a gene involved in efflux of toxic substances (matE) was significantly up-regulated by BBR (1.72 log2FC; P = 0.01), which may indicate export of BBR by this gene product.
BBR treatment significantly affected the topology and complexity of networks between Bilophila gene expression and bile acid profiles in liver, serum, and cecum (Fig. 3B & C). Many of the gene expression networks, sparse in the control, and unconnected with bile acids, become highly interconnected during BBR treatment and relate either directly or indirectly to increased bile acid concentrations. Notably, universal stress protein (uspA) was highly expressed in the BBR group relative to control (3.36 log2FC; P = 0.05) as were genes involved in DNA recombination and repair (sbcC, xseA, uvrD). The expression of uspA revealed a strong positive correlation with cecal bile acids, including DCA (r = 1.0; P = 0.0), TUDCA3S (r = 0.97; P = 0.0), and MDCA (r = 0.97; P = 0.0). DCA also correlated strongly with a recently described glycyl radical enzyme (T370_R50117375; r = 1.0; P = 0.0) involved in sulfide formation from taurine (Fig. 3D) . DCA and the glycyl radical enzyme shared strong positive correlation with NADH dehydrogenase subunit 1 (T370_R50109290; r = 1.0; P = 0.0) as well as direct positive correlations with glycine dehydrogenase (T370_R50113235; r = 1.0; P = 0.0). In control mice, nitrate reductase γ-subunit is positively correlated with cecal CA (r = 0.93; P = 0.001) (Fig. 3B; Supplementary Dataset 1); whereas there are no direct correlations between cecal bile acids and γ-subunit in berberine treatment (Fig. 3C; Supplementary Dataset 1).
Bacteroides uniformis. Seventeen genes were significantly differentially regulated in B. uniformis by BBR (Fig. 4A; Supplementary Dataset 2–4). Two genes were identified whose expression correlated to bile acids: the highly up-regulated NAD(P)H nitroreductase (ERS852554_00867; 2.49 log2FC; P = 0.02), as well as the Na+/H+ antiporter (1.34 log2FC; P = 0.02). A high degree of positive connectivity was observed between total cecal bile acids (r = 0.8; P = 0.046), and primary unconjugated bile acids including β-MCA (r = 0.8; P = 0.046), UCA (r = 0.87; P = 0.015) and 7-oxo-CA (r = 0.8; P = 0.046), as well as expression of acetyl-CoA carboxylase biotin carboxyl carrier protein which was induced by BBR treatment (BLV12_RS03955; 2.47 log2FC P = 4.51E-03; FDR = 0.48) (Fig. 4B & C; Supplementary Dataset 2–4). Chromate transporter (BLV12_RS04500; Log2FC = 3.06; P = 0.01; FDR = 0.58) was negatively correlated with total bile acids in the liver (r = -0.9; P = 0.006) and individual conjugated, sulfated, and primary bile acids (r = -0.9 to -1.0; P = 0.006 to P < 0.001).
Bacteroides vulgatus. BBR differentially altered expression of 105 genes in B. vulgatus (Fig. 5A; Supplementary Dataset 2–4). Most notably, BBR increased the expression of a polycistronic operon encoding predicted multidrug efflux pump subunits—periplasmic adaptor subunit efflux RND (3.72 log2FC; P = 7.23E-08; FDR = 1.23E-05), AcrB/AcrD/AcrF (3.00 log2FC; P = 5.54E-07; 6.87E-05), and TolC (2.71 log2FC; P = 1.6E-07; FDR = 2.26E-05). Expression of the efflux resistance-nodulation-division transporter periplasmic subunit (BVU_RS20445) was positively associated with DCA in the cecum (r = 1.0; P < 0.001) as well as cecal MDCA (r = 0.97; P = 0.0) and TUDCA-3S (r = 0.97; P < 0.001). In addition, a gene encoding a predicted cation/H(+) antiporter was positively correlated with total cecal bile acids (r = 1.0; P < 0.001) (Fig. 5B – D; Supplementary Dataset 1).
BBR induced increased expression of sialidase (1.44 log2FC, P = 2.93E-05, FDR = 1.81E-03) and other genes involved in mucin degradation including N-acetylneuraminate lyase (1.33 log2FC, P = 9.64E-04, FDR = 0.04), N-acylglucosamine 2-epimerase (1.18 log2FC, P = 3.95E-08, FDR = 0.08), α-1,2-mannosidase (0.93 log2FC, P = 1.64E-03 FDR = 0.05), α-L-fucosidase (0.75 log2FC, P = 0.02, FDR = 0.23), and α-1,2-C3/C4-fucosidase (0.70 log2FC; P = 0.02; FDR = 0.23). SusC and SusD outer membrane protein encoding genes, involved in binding and uptake of carbohydrates, were observed in the correlation network in the BBR treated group. In particular, BVU_1844 was differentially expressed in the BBR group (1.18 log2FC; P = 0.03; FDR = 0.31) and negatively correlated with cecal T-β-MCA (r = -0.97; P < 0.001). These data may indicate a “ramping up” of carbohydrate metabolism and may explain the significant increase in total SCFA levels reported previously during BBR intake .
Parabacteroides distasonis. Two operons encoding predicted tryptophan (BDI_RS02910-BDI_RS02940) and leucine biosynthesis (BVU_RS10160-BVU_RS10180; BVU_RS12860-BVU_RS12880) pathways were among the most highly differentially up-regulated genes in P. distasonis in the presence of BBR (Fig. 6A; Supplementary Datasets 2–4). Genes encoding a predicted efflux RND transporter periplasmic adaptor (BDI_RS01695; 2.33 log2FC; P = 1.50E-04; FDR = 0.01) and TolC (BDI_RS00690; 2.30 log2FC; P = 2.04E-06; FDR = 1.37E-03) were also highly expressed, which may reflect adaptation to increased bile salt concentrations in response to BBR. Network analysis showed sparse associations between transcripts and cecal bile acids in control mouse ceca and liver (Fig. 6B & C; Supplementary Dataset 1), but tight interconnections between cecal and liver bile acids and transcripts after BBR treatment (Fig. 6D & E). While TolC expression was not correlated with cecal bile acids, efflux RND transporter had strong positive correlations with cecal DCA (r = 1.0; P < 0.001), TUDCA3S (r = 0.97; P < 0.001), and MDCA (r = 0.97; P < 0.001) in the BBR group (Fig. 6D). Positive correlations were also observed between cecal bile acids and genes involved in leucine and tryptophan biosynthesis (Fig. 6D). Also notable is the positive association between LCA in the cecum and the Na+/H + antiporter NhaA (BDI_R503835; 1.53 log2FC; P = 4.38E-03; FDR = 0.12; r = 0.97; P < 0.001) (Fig. 6D). As in B. vulgatus, several SusC/SusD membrane associated protein genes were also differentially regulated by BBR and correlate directly or indirectly with bile acids in the cecum and liver (Fig. 6D & E). Interestingly, P. distasonis is observed to differentially express multidrug transporter matE (WP_011966429.1; 2.60 log2FC; P = 5.51E-05; FDR = 0.04), which does not correlate with bile acids and may indicate an export protein important for removing intracellular BBR. Indeed, MATE transporters have been shown previously to catalyze xenobiotic compound efflux in a Na + or H + dependent manner [18, 19].
Clostridium hiranonis. The most highly expressed genes in response to BBR in C. hiranonis were a 6 ORF polycistron encoding a helix-turn-helix xenobiotic response protein, chaperone, ATPase, ATPase, MBL fold hydrolase, and dinitrogenase iron-molybdenum cofactor (4.4 to 2.36 log2FC, P = 5.18E-05 to P = 0.024) (Fig. 7A & B; Supplementary Dataset 2–4). Genes involved in peptidoglycan synthesis (murJ and murF) and maintenance of the cell-wall were also significantly up-regulated by BBR. Stress-induced genes, genes involved in DNA repair (recN 0.98 Log2FC; P = 0.021), and the exodeoxyribonuclease VII large subunit (1.19 Log2FC; P = 0.015) were also induced by BBR treatment (Fig. 7A). The expression of exodeoxyribonuclease VII large subunit correlated positively with DCA in the liver (r = 0.95; P < 0.001) and less-positive correlations were observed for other liver bile acids (MDCA, T-β-MCA, TCA, TCDCA) and negatively with 3-dehydro-4,6-CA in the control cecum (r = -0.89; P = 0.008) (Fig. 7C & D; Supplementary Dataset 1). Expression of recN negatively correlated with total cecal bile acids (r = -0.89; P = 0.004) and σ54-dependent Fis family transcriptional regulator (r = -0.9; P = 0.006). Interestingly, treatment with BBR changed this interaction (Fig. 7D; Supplementary Dataset 1). The expression of exodeoxyribonuclease VII large subunit was not contingent on σ54-dependent Fis family transcriptional regulator, but was negatively correlated with serum TCA (r = -0.86; P = 0.017) and total serum bile acids (r = -0.86; P = 0.017). Expression of σ54-dependent Fis family transcriptional regulator correlated positively with liver bile acids (r = 0.8 to 1.0; P values from < 0.05 to < 0.001), and positively correlated with Na+/H+ antiporter which was itself positively correlated with CA in the liver (r = 0.8; P < 0.05) but negatively correlated with numerous cecal bile acids (r = -0.89 to 0.0; P values from 0.007 to < 0.001). It is possible that by importing protonated bile acids, it is not necessary to exchange ions and expression of the Na+/H+ antiporter may decrease.
A putative adhesin was also significantly expressed in the presence of BBR (2.68 Log2FC; P = 2.57E-03). Numerous copies of genes encoding putative cell wall-binding repeat 2 family protein were significantly down-regulated by BBR ranging from − 1.60 log2FC (FDR = 9.04E-04) to -5.59 log2FC (FDR = 7.42E-04). This may indicate modulation of the peptidoglycan layer by BBR treatment. Further support of this was the significant increase in expression of murJ, encoding Lipid II flippase (1.69 log2FC; P = 0.02; FDR = 0.16) and murR transcriptional regulator (1.53 log2FC; P = 1.90E-03; FDR = 0.04) with a trend for murG (1.12 log2FC; P = 0.10; FDR = 0.39) and murF (0.90 log2FC; P = 0.01; FDR = 0.12) was observed. Thus, C. hiranonis gene expression reflects responsiveness to bile acid-induced stress during BBR treatment.
Clostridium hylemonae. BBR differentially regulated 92 genes in C. hylemonae. Of note, a gene predicted to encode the septation ring formation regulator, EzrA, was among the most highly up-regulated genes (2.41 log2FC; P = 1.39E-03) (Fig. 8A; Supplementary Dataset 2–4), but did not correlate with cecal bile acids. Specifically, genes involved in bile acid 7α-dehydroxylation by C. hylemonae were down-regulated regulated, including baiB encoding bile acid coenzyme A ligase (1.45 Log2FC; P = 0.04), and baiCD encoding bile acid NAD-dependent 3-dehydro-4-oxidoreductase (-2.42 Log2FC; P = 3.61E-03) (Fig. 8A). Phage genes, including holin (2.10 log2FC; P = 2.8E-04) and siphovirus DUF859 (1.96 Log2FC; P = 1.27E-3) as well as type I-C CRISPR Cas8c/Csd1 (1.09 log2FC; P = 0.04) were up-regulated by BBR. In control mice, phage holin expression was positively correlated with cecal DCA (r = 0.81; P = 0.021), but negatively correlated with cecal T-β-MCA (r = -0.83; P = 0.016) and β-MCA in the liver (r = -0.94; P = 0.0) (Fig. 8B; Supplementary Dataset 1). In BBR treated mice, phage holin was positively associated with total cecal bile acids (r = 1.0; P = 0.0) and had a strengthened positive correlation with cecal DCA (r = 0.9; P = 0.006) (Fig. 8C & D; Supplementary Dataset 1).
Cecal RNA-Seq analysis revealed two polycistronic operons involved in the Stickland fermentation of glycine, including the glycine dehydrogenase and glycine reductase pathway genes and the formation of cofactors such as lipoate that were significantly up-regulated by BBR (Fig. 8A). In control ceca, expressed genes involved in glycine reductase (FolD, grdD) appeared to be indirectly and negatively correlated to bile acids via transcription terminator/antiterminator NusG (Fig. 8B). Lipoate-protein ligase A expression was positively correlated with expression of metabolic genes as well as CA-4,6-3-one. In BBR treated mice, lipoate-protein A displayed strong positive correlations with total cecal bile acids (r = 0.8; P = 0.046) and individual cecal bile acids, such as DCA (r = 0.9; P = 0.006) as well as weak negative correlations with liver bile acids (Fig. 8C). Positive correlations were observed between lipoate-protein ligase A and FolD, glycine cleavage protein T, and dihydrolipoyl dehydrogenase, indicating that the significant increase in glycine metabolism with BBR treatment was at least partially driven by increased bile acid concentration in the cecum.
Genes involved in sporulation in C. hylemonae including spore coat associated protein (cotJA; 2.29 log2FC; P = 2.89E-04), N-acetylmuramoyl-L-alanine amidase (cwlD; 1.54 log2FC; P = 0.03), spore germination protein (1.43 log2FC; P = 0.03), acid-soluble spore protein (0.96 log2FC; P = 0.04) were observed in response to BBR treatment. Acid-soluble spore protein was negatively correlated in control mice with liver bile acids, whereas BBR treatment resulted in a positive correlation with liver bile acids. This may indicate that up-regulation of genes involved in cell-wall maintenance and metabolism reflects the effects of both BBR and bile acids.