Taxonomic profiling and populational patterns of bacterial bile salt hydrolase (BSH) genes based on worldwide human gut microbiome

The molecular weights of the BSH subunits ranged from 28 to 50 kDa [30]. Comparison of sequences revealed that the number of amino acids in BSHs from Clostridium, Bifidobacterium, Lactobacillus, and Enterococcus were 329, 316, 325, and 326, respectively (Additional file 1: Figure S4a). Pairwise amino acid sequence alignments of the BSHs showed that the highest identity was 54.29% (comparison of BSHs from Lactobacillus and Enterococcus), and the lowest identity was 34.91% (comparison of BSHs from Bifidobacterium and Enterococcus) (Additional file 1: Figure S4b). Overall, their average value of sequence identity was 40.41% (Additional file 1: Figure S4a).

The secondary structures (the pattern of hydrogen bonds between the amino hydrogen and carboxyl oxygen atoms in the peptide backbone) of these BSH showed that the six reported active site residues were conservative (Fig. 1a, highlight by red frame), although the number of α-helix and β-strand were different. The topological diagram of BSH showed that the domain had a six-layered structure of composition βαββαβ, and the core of BSH was composed of two sandwiched antiparallel α-sheets, which were conformed to the structure characteristics of the Ntn-hydrolase family [29, 31] (Fig. 1b).

According to previous studies, the conserved and functional important residues in these four BSH were Cys2, Arg18/18/16/16, Asp21/21/19/19, Asn82/82/79/79, Asn175/173/170/171, and Arg228/226/223/224 [26–29]. These six residues were closed in the 3D structure to form the core active site of BSH although they were far apart in BSH (Fig. 1c). It is noteworthy that the functional active sites in BSHs are not only these six residues. For instance, there were totally 30 residues that are contacted with various substrates in BSH of 2BJF (Additional file 1: Figure S5, Additional file 2: Table S1).

These results indicated that although there exist much sequence differences among different BSH, the functional residues, αββα fold of secondary structure, and core active site are conserved.

The homology of BSHs between genera demonstrated that they were mainly distributed in two intervals, one at about 35% and the other at about 50% (Fig. 2e). However, there was also a lower identity of BSHs within the genus at about 50% (Fig. 2e). These results indicated that distinguishing BSHs by genera might not be a rational method because of the paralogs of BSHs harbored in the bacterial genome.

The BSHs obtained from the HMP database were further screened in the gut microbiome of worldwide populations (Additional file: Table S3). Among the 561 BSHs we obtained from the HMP database, only 156 were presented in the gut microbiome of 591 people from 11 different countries of six continents (Additional file: Table S4). The drastic decrease in the number of strains (591 to 156, 26.40%) indicates that many BSH-expressing bacteria are very low in abundance in the human gut microbiome and undetectable under the current sequencing depth. 2 2

Based on the cumulative RA distribution of BSHs in different populations (Fig. 3d, Additional file 2: Table S4), the America (US) population exhibited the highest cumulative RA of BSHs (3.74 × 10⁻⁴), while the Hadza ethnic group of Tanzania (HZ) exhibited the lowest (6.20 × 10⁻⁵). The lower cumulative RA of BSHs in the HZ population may have been because they live in habitats almost completely isolated from other humans and have had relatively little modification to their basic way of life for hundreds of years [33]. Moreover, majority of BSHs were identified from five genera, Bacteroides, Blautia, Eubacterium, Clostridium, and Roseburia, which represented about 71.31% of the total abundance of BSHs in human gut metagenomes (Fig. 3d, Additional file 2: Table S5).

Multivariable regression analysis was implemented to evaluate the relationship of BSHs abundance with gender, age, BMI, and population (Fig. 3). Gender, age, and BMI were not significantly correlated with the abundance of BSHs after adjustment (Fig. 3a, b, and c). However, significant associations were found between the cumulative RA of BSHs and populations from different geographical regions (p^b = 1.2e-28, Fig. 3d). The populational factor explained the highest variations (35.15%) of BSHs abundance among the four factors (Fig. 3e). These results were consistent with the observation that disparity of the gut microbiota composition in different populations was mainly caused by geography [34].

Among the genera distribution (middle panel of Fig 4, Additional file 2: Table S4) and population patterns (right panel of Fig. 4, Additional file 2: Table S4) of each phylotype, BSH-T0 contained only 0.25% of the RA of BSHs, which were not shown in the gut microbiome of HZ population and distributed in Clostridium, Intestinibacter, Lactobacillus, and Enterococcus, respectively. BSH-T1 contained 38.04% of the total RA of BSHs, which were mainly distributed in Eubacterium, Blautia, Clostridium, Roseburia, and Ruminococcus and could be found in the gut microbiota of all 11 populations. BSH-T2 contained 1.01% of the total RA of BSHs, and the major genera of this phylotype were Streptococcus and Enterococcus. Notably, this phylotype was not found in gut microbiota of the population of China (CN; right panel of Fig 4, Additional file 2: Table S6). Sequences of BSH-T3 were all from Lactobacillus, which were only found in CN, Japan (JP), Austria (AT), and France (FR) populations and showed higher RA in CN and AT populations (Fig. 4, Additional file 2: Table S6). BSH-T4 represented 2.74% of the total RA of BSHs, which were mainly distributed in Bifidobacterium and Collinsella and could not be found in HZ population (Fig. 4, Additional file 2: Table S6). The RA of BSH-T4 in gut microbiota of the population of JP was higher than that of other populations (Fig. 4, Additional file 2: Table S6). BSH-T5 contained 23.63% of the total RA of BSHs, which were mainly distributed in Bacteroides (Fig. 4). BSH-T6 contained 32.2% of the total RA of BSHs, which were also mainly distributed in Bacteroides (Fig. 4). However, there were no BSHs from this phylotype could be found in HZ and Peru (PE) population (Fig. 4, Additional file 2: Table S6). BSH-T7, which contained 1.98% of the total RA of BSHs mainly distributed in Blautia and could be found in the population from all countries (Fig. 4, Additional file 2: Table S6).

There were 28 strains with BSHs paralogs in 120 strains, the BSHs from Bacteroides, Clostridium, Lactobacillus, Ruminococcus, and Marvinbryantia showed varied phylotype distributions. Overall, nine strains with paralogs were distributed in the same phylotype (Additional file 1: Figure S8, marked by black triangles), 15 strains were distributed in two phylotypes (Additional file 1: Figure S8, marked by red triangles). Specifically, most of the strain with paralogs were from Bacteroides, and mainly distributed in two phylotypes, i.e., BSH-T5 and BSH-T6 (right panel of Fig. 4, Additional file 1: Figure S8). Four strains were distributed both within and between phylotypes (Additional file 1: Figure S8). Thus, different paralogs of BSHs within a genus could have sequence dissimilarity, which might lead to variable functional roles of these two genera in bile acid metabolism.

Representative BSHs of the eight phylotypes were synthesized and purified (Additional file 1: Figure S9a) as described in the “Methods” section. The deconjugation of BSH at lowest concentration (0.1 mM) of substrates was confirmed by LC-MS/MS (right panel of Fig. 5, Additional file 2: Table S9). The results of enzyme activity assay shown that BSH-T0, BSH-T1, BSH-T3, BSH-T4, and BSH-T7 performed higher specific activity (the activity of an enzyme per milligram of total protein, using the highest enzyme activity as 100% for each substrate) among seven bile salts (middle and right panel of Fig. 5, Additional file 2: Table S8). However, different BSH phylotypes display selective deconjugation activity based on substrates. The specific activity of BSH-T5 was lower when the substrate was glycocholic acid (GCA) and both GCA and taurocholic acid (TCA) for BSH-T6 (middle and right panel of Fig. 5, Additional file 1: Figure S18). In particular, BSH-T1, BSH-T3, and BSH-T4 showed the highest specific activity with GCA, while BSH-T0 and BSH-T2 showed the highest specific activity with glycochenodeoxycholic acid (GCDCA), and BSH-T0 showed the highest specific activity when taurochenodeoxycholic acid (TCDCA) was the substrate (Additional file 1: Figure S18). It is worth noting that BSH-T1, which had the highest abundance of BSHs in the human gut microbiota, exhibited higher deconjugation activity (middle and right panel of Fig. 5). BSH-T2 showed lower enzyme activity both in silico and in vitro (Fig. 5, Additional file 1: Figure S12, S18). Comparatively, BSH-T7 showed higher deconjugation activity in vitro but not in silico (middle and right panel of Fig. 5, Additional file 1: Figure S17, S18). This discrepancy was likely because the computational molecular docking only partially reflected the actual enzyme activities. Nevertheless, computational work is helpful to understand the biological activity at the molecular level.

To further investigate the functional implications of BSH, we analyzed the relationship between the RA of BSHs and disease status of different groups, including populations of geographical regions and health/disease cohorts.

First, we analyzed the relationship between the RA of BSHs in our target populations with the World Health Organization (WHO) released phenotypes, namely, death rate of diabetes, death rate of cardiovascular diseases (CVD), mean blood cholesterol, and BMI of obesity (Additional file 1: Figure S19, Additional file 2: Table S10). We found that the RA of all BSH was significantly correlated with death rate of diabetes (r = − 0.65, p = 0.03) (Additional file 1: Figure S19a, entry 7 of Additional file 1: Figure S19c), and the RA of all BSH-T0 was significantly correlated with death rate of CVD (r = 0.92, p = 0.0041) (Additional file 1: Figure S19b, entry 7 of Additional file 1: Figure S19c). These findings might indicate that the RA of BHSs was relevant to CVD risk among different populations. However, it should be noted that the metagenome data and epidemiological data did not originate from the same cohort in this association study, which may be a limitation of this investigation. Nevertheless, these results still shed light on the relationship of BSH with various diseases.

The known 3D structures of BSHs were searched in the protein data bank (PDB) [25] using the keyword “bile salt hydrolase” and “Choloylglycine hydrolase”. Only 12 of the 25 searched results are our purpose structures (Additional file 1: Figure S2b). Pairwise amino acid sequence alignments of the BSHs were performed by BLASTP (v 2.2.29+) [64], and multiple alignments were conducted using DNAMAN (v 8.0). The secondary structures of these BSHs were downloaded from PDBsum [65]. Comparison of 3D structures was conducted using MOE (v 2014).

The reference genomes of 1751 bacterial strains covering 1253 species were obtained from the HMP database [32] in September 2014. The genes and related proteins from these bacterial genomes were predicted by MetaGeneMark (v 2.8) [66], and the taxonomic information regarding these genes and proteins was directly extracted from the strain names.

The query BSH protein sequences were collected from the Refseq database of NCBI database [24] using the keyword “bile salt hydrolase” and “Choloylglycine hydrolase”. Then, the total number of amino acids in the sequences was limited between 300 and 400. The sequence identity interval of the BSHs between and within genera was calculated from these 765 sequences to define the thresholds. Thus, the final 591 BSHs reference sequences (screening was performed by controlling the total length of the sequences from 300 to 400) were identified from HMP database by taking the initial 765 BSH as query and using BLASTP with e-value of 1e-5 and sequence identity of 45% as cutoff (Additional file 1: Figure S2a).

The metagenomic sequence data of individuals were collected from 11 populations of six continents, including the Hadza ethnic group of Tanzania (HZ) of Africa [33, 67], China (CN) of Asia [50, 68], Japan (JP) of Asia [69], South Korea (KR) of Asia [70], Denmark (DK) of Europe [71], Sweden (SE) of Europe [72], Austria (AT) of Europe [56], France (FR) of Europe [73], Australia (AU, PRJEB6092) of Oceania, the United States (US) of North America (N-America) [21, 32], and Peru (PE) of South America (S-America) [74].

To construct metagenomic datasets of healthy individuals from each country, we screened out the data for individuals with 18.5 < BMI < 29.9 kg/m² [69, 75], 12 < age < 75 [69, 76], and those designated with diseases that were excluded from the data. Although we could not access the metadata for individuals from the United States, we used all data from healthy individuals with an average BMI of 24 ± 4 kg/m² for this cohort. Finally, a total of 581 healthy individuals were selected from these 11 countries for analysis (Additional file 2: Table S3).

To study the differences in the BSH relative abundance (RA) between healthy individuals and patients, we used metagenomic datasets of individuals from different countries who were identified with diseases expected to be related to deviations in bile acids such as colorectal cancer (CRC), adenoma, type 2 diabetes (T2D), impaired glucose tolerance (IGT), and atherosclerosis (AS) [58] (Additional file 2: Table S3).

All raw sequencing reads were assessed and filtered using the FASTX-Toolkit, and the high quality microbiome sequencing reads were assembled with the SOAPdenovo2 (v 2.04) [77] package. After assembly, contigs with at least 500 bp were further used to predict the genes with MetaGeneMark (v 2.8), after which a non-redundant protein set was constructed by pair-wise comparison of all protein sequences within populations using BLASTP (v 35x1) [78] with 95% identity and 90% overlapping thresholds. The relative abundance (RA) of each protein sequence in each individual was calculated based on the number of read pairs mapped to the gene over the length of the protein sequence and divided by the summary of sequence abundance per individual [50]. The cumulative RA was calculated by the sum of RA of each genus in each population. The BSHs were identified in the above population database using BLASTP with an e-value = 1e-5 and 62% sequence identity as the cutoff values.

A phylogenetic tree was built using the maximum likelihood method in the MEGA software (v 7.0). Dendroscope (v 3.4.7) was used to embellish the phylogenetic tree by adjusting the labels and filling the colors as needed.

Nationwide blood cholesterol (2008), obese BMI (2014), diabetes death rate (per 100,000 individuals; 2012) and cardiovascular diseases death rate (per 100,000 individuals; 2012) data by country were downloaded from the World Health Organization (WHO) [79] (Additional file 2: Table S10).

The online software, Protein Homology/analogY Recognition Engine (Phyre2, V 2.04) [80], was employed to predict the homologous structure of BSHs using intensive mode [81]. Ligands including glycocholic acid (GCA, C₂₆H₄₃NO₆), glycochenodeoxycholic acid (GCDCA, C₂₆H₄₃NO₅), glycine deoxycholic acid (GDCA, C₂₆H₄₃NO₅), taurocholic acid (TCA, C₂₆H₄₅NO₇S), taurochenodeoxycholic acid (TCDCA, C₂₆H₄₅NO₆S), taurodeoxycholic acid (TDCA, C₂₆H₄₅NO₆S), and tauroursodeoxycholic acid (TUDCA, C₂₆H₄₅NO₆S) were obtained from the ZINC database [82]. The AutoDock program (v 4.2.6) [83] was employed to generate an ensemble of docked conformations for each ligand bound to its target. Four units of BSH were found to bond with single bile acid substrate molecules by forming a tetramer. We then used the genetic algorithm for conformation searches and conducted 100 individual genetic algorithm (GA) runs to generate 100 docked conformations for each ligand.

Sodium glycocholate hydrate (CAS: 863-57-0), sodium glycochenodeoxycholate (CAS: 16564-43-5), sodium glycodeoxycholate (CAS: 16409-34-0), sodium taurocholate hydrate (CAS: 345909-26-4), sodium taurochnodeoxycholate (CAS: 6009-98-9), sodium taurodeoxycholate hydrate (CAS: 207737-97-1), sodium tauroursodeoxycholate (CAS: 35807-85-3), CA-d4 (CAS: 116380-66-6), DCA-d4 (CAS: 112076-61-6), CDCA-d4 (CAS: 99102-69-9), UDCA-d4 (CAS: 347841-46-7), glycine (CAS: 56-40-6), and taurine (CAS: 107-35-7) were used in the enzyme assay. pET28a (+) was used for the expression of His-tagged (6x) recombinant BSH (Genscript, Nanjing, China). The target BSH sequences were ligated into the NcoI/XhoI-digested vector pET28a (+) (T7 promoter), resulting in six recombinant plasmid pET28a (+)-bsh genes being transformed into E. coli BL21 (DE3) competent cells.

The eight bsh-recombinant bacteria were inoculated into Luria-Bertani (LB) broth containing 50 μg·ml⁻¹ kanamycin, and expression of the BSH genes was induced by the addition of isopropyl β-D-thiogalactoside (IPTG, 0.5 mM). The cell pellet obtained by harvesting was subsequently resuspended (1:10) in buffer solution (0.05 M Tris, 0.05 M NaCl, 0.5 mM EDTA, 5% Glycerol, pH 7.9) and disrupted by sonication (alternating pulses: on for 3 s, off for 3 s; 60% amplitude) until clear solutions were obtained.

The BSH was purified by nickel-nitrilotriacetic acid (Ni²⁺-NTA) agarose column (HisTrap™ HP, GE Healthcare, USA). The presence and purity of BSH in each sample were confirmed by SDS-PAGE. After which pure BSH proteins were stored at − 80 °C after subsequently lyophilized.

The BSH specific activity was determined by measuring the release of amino acids (glycine and taurine) from conjugated bile salts. The amounts of amino acids liberated by BSH were determined by ninhydrin assay [84]. Protein concentrations were determined using a Bradford Protein Assay kit (Beyotime, Nanjing, China). The purified BSH-lyophilized powders were diluted with 20 mM phosphate buffer (pH 6.5) until the concentration of BSH was 0.65 mg·ml⁻¹. Next, 10 μl of BSH liquid, 10 μl of bile salts, and 180 μl of reaction buffer (20 mM phosphate buffer, pH 6.5) were mixed, and 10 μl of liquid paraffin were added. The samples were subsequently incubated at 37 °C for 30 min, after which used 200 μl of 15% (w/v) trichloroacetic acid to terminate the reaction, and the sample was centrifuged to remove the precipitates. Mix 10 μl of reaction supernatant with 190 μl of ninhydrin reagent and kept in a boiling water bath for 15 min. The absorbance at 570 nm was measured in a 96-well plate. A standard BSH activity curve was subsequently prepared for glycine and taurine (Additional file 1: Figure S11b). One unit of BSH activity was defined as the amount of enzyme that released 1 μmol of amino acid from the substrate per min.

Analysis was performed on a triple quadrupole tandem liquid chromatography-mass spectrometry (LC-MS/MS) system (LCMS8050, Shimadzu). Chromatographic separation was performed on an ACQUITY UPLC BEH C18 (2.1 × 100 mm, 1.7 μm) column. The mobile phase consisted of 0.1% formic acid (mobile phase A) in water and 0.1% formic acid in acetonitrile (mobile phase B) running at a flow rate of 0.4 ml/min. The gradient elution program was 5% B at 0–1 min, 5–42% at 1–2.5 min, 42–45% at 2.5–5.5 min, 45–60% at 5.5–8 min, 60–95% at 8-9 min, 95% at 9–9.5 min, back to initial conditions, and 2 min for equilibration. The column was maintained at 55 °C and the injection volume of each sample was 1 μl.

The LC-MS/MS system control and data analyses were performed by LabSolutions software (the software version: version 5.65). The ion source parameters were set as follows: nebulizing gas flow of 2 l/min, heating gas flow of 10 l/min, interface temperature of 300 °C, DL temperature of 250 °C, heat block temperature of 400 °C, and drying flow of 10 l/min. The data were collected with multiple reaction monitor (MRM) in negative mode.

An isotope-labeled standard calibration approach with standard addition was used to avoid the matrix effects and ensure the accuracy of measurement. Calibration curves were constructed each day using seven calibrators prepared from pooled plasma spiked with CA, DCA, CDCA, and UDCA at concentration range of 0.00625–0.2 mM and CA-d4, DCA-d4, CDCA-d4, and UDCA-d4 at 0.003125 mM (Additional file: Figure S21, S22). 1

All values were expressed as mean ± SEM. Statistically significant differences between two groups were determined by Mann-Whitney U test followed by false discovery rate correction. The relation analysis was done with the Spearman rank correlations. Multiple comparisons were performed by multivariable-adjusted analysis using linear regression model. All analyses were performed using R (v 3.3.2), and p < 0.05 was considered as statistically significant.

Additional file 1: Figure S1. The enterohepatic circulation of bile acids. Figure S2. The analysis flowcharts. Figure S3. Genera characterization of BSHs obtained from NCBI database. Figure S4. The sequence identity of reported BSH sequences. Figure S5. Reported conserved residues of BSH in PDB database. Figure S6. Taxonomic identification of BSHs in HMP database. Figure S7. The Taxonomic identification and paralogs distribution of BSHs in HMP. Figure S8. Phylogenetic tree of BSHs in the gut microbiome of 11 populations. Figure S9. SDS-PAGE of BSH and Standard curve used in ninhydrin assay. Figure S10. Molecular docking results between BSH-T0 and bile acids. Figure S11. Molecular docking studies between BSH-T1 and bile acids. Figure S12. Molecular docking studies between BSH-T2 and bile acids. Figure S13. Molecular docking studies between BSH-T3 and bile acids. Figure S14. Molecular docking studies between BSH-T4 and bile acids. Figure S15. Molecular docking studies between BSH-T5 and bile acids. Figure S16. Molecular docking studies between BSH-T6 and bile acids. Figure S17. Molecular docking studies between BSH-T7 and bile acids. Figure S18. Enzyme activity comparisons between substrates in each BSH phylotype. Figure S19. Relationship between average RA of BSHs with WHO released phenotypes in 11 populations. Figure S20. The relative abundance of BSHs among case-control cohorts. Figure S21. Standard curve of bile acids used in LC-MS/MS. Figure S22. Typical LC spectrum of different bile acids. (PDF 4413 kb) Additional file 2: Table S1. Reported conserved residues of BSH in PDB database. Table S2. Taxonomic identification of BSHs in HMP database. Table S3. The information of individuals from the 11 populations. Table S4. Cumulative relative abundance of reclassified BSHs in 11 populations. Table S5. Cumulative relative abundance of BSHs in different genera in 11 populations. Table S6. The percentage of relative abundance of each BSH phylotype and each population. Table S7. Information of represent BSHs used in molecular docking and enzyme activity assay. Table S8. BSH activity at a substrate concentration of 20 mM by ninhydrin assay. Table S9. Percent deconjugation of BSH at a substrate concentration of 0.1 mM. Table S10. Information of WHO released phenotypes in 11 populations. (XLSX 100 kb)