Harnessing a germ‑free mouse gut bioreactor for directed evolution of probiotics to combat non-alcoholic fatty liver disease

Domain: chemistry
Published: 2026-02-24
Editorial outcome: major_revision

Abstract
Directed evolution has transformed biomolecular engineering but remains largely untapped for probiotic optimization, despite its immense promise for human health maintenance and disease therapy. Here, we present an in vivo, host-mediated directed-evolution platform that harnesses the gut’s endogenous selective pressures to drive functional enhancement of probiotics. Using Bifidobacterium animalis subsp. lactis as a model, we expose germ-free male mice to stepwise increases in bile-acid stress via a high-fat, high-cholesterol diet. Compared to in vitro evolution, which fails to produce any functional gains, our host-mediated approach yields a variant exhibiting a 77% increase in bile acid metabolism. Multi-omics analysis identifies two critical single-nucleotide variants (SNVs) simultaneously: one in the upstream region of cbh, encoding bile salt hydrolase, and a non-synonymous mutation in mdr, a bile-acid efflux transporter. Functional validation assays confirm that these mutations drive increased corresponding gene expression and enhance substrate binding efficiency. Moreover, to demonstrate its translational potential, we administer the parental and adapted strains daily in a male diet-induced mouse model of non-alcoholic fatty liver disease (NAFLD). We find that the adapted strain significantly improves bile-acid homeostasis, reduces hepatic steatosis, lowers inflammatory and lipid biomarkers, and enhances histological liver health compared to the parental strain. Our work establishes the host gut as a living evolutionary bioreactor for precision engineering of probiotics, enabling targeted phenotypic enhancement in vivo through natural selection. This framework paves the way for personalized, functionally tailored microbiome therapeutics and sets a foundation for next-generation live biotherapeutic agents.

Introduction
Directed evolution is an iterative process of mutagenesis and phenotypic screening that enables efficient optimization of biomolecular functions for specific traits. It harnesses evolutionary principles by introducing genetic diversity under selective pressure, driving the emergence of elite mutants through progressive screening cycles1. Described by the Nobel Committee as harnessing evolution for the greatest benefit to humankind, this methodology has found wide transformative applications in enzyme engineering, antibody development, and metabolic pathway refinement2–4. Currently, the most widely used strategies for generating targeted microbial adaptation include two major strategies: (1) induction-based methods, such as adaptive laboratory evolution (ALE)5–7; and (2) gene-editing-based methods, such as CRISPR-Cas9-based genome editing7,8. For example, short-term ALE using phenolic acids (such as ferulic acid, p-hydroxybenzoic acid, and vanillic acid) as model inhibitors yielded Saccharomyces cerevisiae strains with 2- to 3-fold increased growth rates and enhanced tolerance to lignin-derived phenolics9. Similarly, as an example of gene editing, CRISPR-Cas9-mediated modification of repetitive elements in the phytosterol biosynthesis pathway of Solanaceae plants enabled a fivefold increase in 7-dehydrocholesterol (7-DHC), a vitamin D₃ precursor, in leaves and immature fruits without impairing plant development10.

Probiotics have garnered considerable attention in recent years due to their multifaceted roles in maintaining gut homeostasis11, regulating metabolic functions12, and reinforcing immunological barriers13. Nevertheless, most studies to date concentrate on screening native strains14 and characterizing their phenotypic functions15 or probiotic properties16, while systematic optimization and directed evolution of their functional potential remain largely unexplored. A few reports have applied ALE in vitro, exposing strains to prolonged passaging under stressors such as acidity17, osmotic pressure18, or antibiotics19 to select for mutants with improved survival or functional traits. Likewise, gene-editing approaches (e.g., CRISPR-Cas) have been used to reconfigure metabolic pathways to enhance gut colonization20, target specific functions21, or boost metabolite production22. However, these in vitro strategies rely on low-throughput screening pipelines that limit scalability in probiotic engineering23. Moreover, strains evolved under laboratory conditions frequently fail to achieve ecological adaptation and functional remodeling when reintroduced into authentic host-gut environments.

The host intestinal environment serves as a natural evolutionary arena that drives the accumulation of adaptive mutations in gut bacteria24,25. Studies demonstrate that gut microbiota acquire single-nucleotide variants (SNVs) at rates orders of magnitude higher in vivo than in vitro (reaching billions in scale) due to intestinal selection pressures26. For instance, genomic analysis of Lactiplantibacillus plantarum HNU082 isolates after 28-day intestinal colonization in model organisms (human, zebrafish, and mouse) revealed a significant enrichment of SNPs associated with acid resistance, membrane protein synthesis, and carbon source metabolism. In contrast, serially passaged in vitro counterparts exhibited no detectable mutations25. This divergence underscores the unique evolutionary pressures exerted by host intestinal microenvironments. Likewise, a parallel comparison showed that Lactiplantibacillus plantarum P-8 passaged in vitro (glucose-limited or standard media) exhibited no genetic changes, while isolates recovered from human feces (in vivo) carried multiple SNPs, including mutations in mucus-adhesion protein and pyruvate oxidase genes27. Although these findings demonstrate that host-imposed intestinal selection pressures can drive functional optimization of probiotics, most studies remain observational and lack phenotype-driven directed-evolution strategies. We therefore seek to investigate whether the host gut itself could serve as a bioreactor for the directed evolution of probiotic strains.

To validate the efficacy of leveraging the host intestinal microenvironment as a bioreactor for directed evolution, we first select Bifidobacterium animalis subsp. lactis H22B656, a probiotic strain with intrinsic bile acid-metabolizing capacity, as our model. We then implement two induction strategies: (i) an in vivo approach, in which germ-free mice are mono-colonized with H22B656 and maintained on a high-fat, high-cholesterol diet (studies have shown that dietary fat intake correlates positively with bile acid secretion and expands the bile acid pool, thus imposing increased bile acid-metabolic pressure on the gut28,29); and (ii) an in vitro approach, consisting of serial passaging in conventional liquid medium supplemented with bile acids. Although the in vitro strategy fails to produce any beneficial phenotypic shift, the in vivo, host-mediated directed evolution yields a strain with a 77% increase in bile acid-metabolizing capacity. Subsequently, we use whole-genome sequencing (WGS) and transcriptomic profiling to reconstruct the evolutionary trajectory of the adapted strain, pinpointing SNVs underlying enhanced bile acid metabolism. We functionally validate two key SNVs (rs1846941 and rs402499) to demonstrate their causal role in modulating bile acid-metabolizing capacity, thereby elucidating the molecular mechanisms of in vivo adaptation. Finally, to evaluate translational potential, we employ a diet-induced nonalcoholic fatty liver disease (NAFLD) mouse model and administer the adapted strain as a dietary intervention, demonstrating its metabolic benefits and affirming its biosafety in vivo. Our study not only introduces an in vivo, host-mediated directed evolution strategy for probiotic optimization but also establishes an approach for targeted interventions against metabolic diseases such as NAFLD, laying a critical foundation for the future clinical translation of disease-tailored probiotic therapies.

Results
Screening of native strains with bile acid-metabolizing capacity
To evaluate the efficacy of the in vitro and in vivo induction strategies for directed probiotic adaptation, we first sought to identify a strain with bile acid-metabolizing capacity as a model. A panel of 12 B. animalis subsp. lactis strains were randomly selected from our library and systematically assessed across four dimensions: (1) a calcium precipitation assay, (2) specific enzyme activity quantification, (3) stress tolerance testing, and (4) growth curve analysis (Fig. 1a; “Methods”). Specifically, to assess bile salt hydrolase (BSH) activity, we plated strains on bile acid-supplemented agar and measured the diameter of the calcium precipitation zone after 48 h. Strains H22B484, H22B525, and H22B656 produced the largest zones (averaging approximately 2.7 cm), indicating robust enzymatic activity (Fig. 1b and Supplementary Fig. 1a). Next, we quantified BSH-specific activity using a ninhydrin-based colorimetric assay. H22B656 exhibited the highest activity (0.4561 U/mg), followed by H22B484 (0.3570 U/mg) and H22B493 (0.3392 U/mg) (Fig. 1c). For stress tolerance assessment, candidate strains were evaluated under both acidic and bile acid conditions to reflect gut survival potential. Under acidic conditions (pH 2.0 for 2 h), seven strains (H22B748, H22B938, H22B801, H22B525, H22B391, H22B648, and H22B656) maintained viability above 25%. When exposed to bile acids (0.3% w/v for 4 h), H22B656 again exhibited the highest survival rate (83.49%), followed by H22B723 (69.50%) and H22B525 (64.67%) (Fig. 1d). Finally, we measured growth kinetics by tracking OD₆₀₀ every 4 h over 36 h. H22B391 demonstrated the highest biomass accumulation (OD₆₀₀=0.5521 at 24 h), followed by H22B656 (0.5031) and H22B748 (0.4778) (Fig. 1e).

To comprehensively evaluate probiotic potential in bile acid metabolism, we integrated these five key indicators from four dimensions and conducted a percentile-based ranking using a multi-dimensional linear summation score for each strain. The resulting composite scores were visualized as radar plots (Supplementary Fig. 1b). Among the candidates, H22B656 achieved the highest composite score as evidenced by the largest radar plot area (425), followed by H22B525 (375) and H22B391 (366.6667) (Supplementary Fig. 1c). Based on this integrative assessment, H22B656 was selected for subsequent adaptive evolution experiments.

The in vivo induction markedly enhanced the bile acid-metabolizing capacity
To enhance the bile acid-metabolizing capacity of B. animalis subsp. lactis H22B656, we employed both in vitro and in vivo induction strategies and systematically evaluated their efficacy in driving targeted adaptive evolution. Using an in vitro induction protocol, H22B656 was serially passaged under 0.3% (w/v) conjugated bile acid pressure for ten days, with transfers every 24 h into fresh, modified MRS-cys medium containing the same bile acid concentration (“Methods”). On day 10, we isolated 60 clones from the induced population and measured their BSH-specific activities (Fig. 2a). Under 0.3% bile acid stress, none of the induced clones exhibited appreciable increases in BSH activity compared to the original strain, with the average specific activity of induced clones remaining ≤1-fold relative to the original (Fig. 2b).

The limited effectiveness of the in vitro induction likely arises from its inability to replicate the complex selective pressures of the gut, particularly the immunological interactions between host and microbiota that drive adaptive evolution25. To address this, we explored an in vivo, host-mediated induction strategy to harness selective forces at the host-microbiota interface and promote specific microbial functions (“Methods”). Using germ-free C57BL/6J mice as a bioreactor to eliminate interference from resident microbes, we imposed a stepwise high-fat diet to progressively intensify bile acid metabolic stress. Over seven weeks, B. animalis subsp. lactis H22B656 colonized mice (200 µL of 10⁹ CFU administered by oral gavage once daily for three consecutive days before diet intervention) received: (1) 30% fat + 0.2% cholesterol for weeks 1–3; (2) 45% fat + 0.2% cholesterol for weeks 3–5; and (3) 60% fat + 0.2% cholesterol for weeks 5–7 (Fig. 2c). To assess the in vivo induction’s efficacy, we also isolated 60 clones from fecal samples collected at different time points and measured their BSH specific activity. Fifteen clones exhibited a significant increase in BSH activity compared to the original (fold change > 1.35), yielding a directed induction success rate of 25% (Fig. 2d). Among these, strain W5S9 exhibited the highest specific activity (fold change = 1.77).

Compared with the in vitro approach, the in vivo strategy imposed a substantially greater selective impact on bile acid metabolism, as reflected by a markedly higher variance in BSH activity (SD = 0.3440 in vivo versus 0.0784 in vitro). This suggests that host-derived pressures can dramatically reprogram microbial metabolic traits (both positively and negatively), generating a diverse pool of adapted candidates. These results confirm the feasibility and utility of in vivo induction as a powerful approach to optimize probiotics’ beneficial effects.

Genetic adaptations in bile acid-metabolism via in vivo induction
We hypothesized that the enhanced bile acid-metabolizing capacity of strain W5S9 resulted from the accumulation of genetic variants that altered key corresponding metabolic pathways. To understand these molecular mechanisms, we conducted whole-genome resequencing (WGS) and transcriptomic profiling (RNA-seq) of W5S9. These analyses aimed to systematically map the single-nucleotide variants (SNVs) acquired during in vivo induction and to assess their impact on gene expression programs governing bile acid metabolism (“Methods”).

To improve the accuracy of mutation identification, we applied two independent computational pipelines (inStrain30 and BCFtools31) and used the intersection of their outputs to define a high-confidence SNV set. This approach identified 282 SNVs, which were then annotated via inStrain (Fig. 3a; Supplementary Fig. 1d; Supplementary Data 1; “Methods”). Based on genomic location, we divided these SNVs into coding and noncoding categories. Of the 244 coding-region SNVs, 39 were synonymous (affecting 38 genes), and 205 were nonsynonymous (spread across 164 genes). As for the remaining 38 SNVs that fell within noncoding regions, promoter-prediction analysis indicated that 11 lie in putative regulatory elements, potentially influencing transcription of nine downstream genes (“Methods”). In addition, given the close genomic linkage, we suspect that these 38 noncoding SNVs may collectively affect the expression of up to 32 downstream genes32 (Fig. 3b).

To explore the phenotypic and functional impact of these mutations, we conducted differential transcriptomic analysis between three adapted clones and three original strains (FDR < 0.001 and |log₂FC| ≥ 0.4 were defined as significant). Among the 164 protein-coding genes harboring non-synonymous SNVs, 16 were significantly upregulated, and 32 were significantly downregulated (Fig. 3c). Of the 41 putative target genes associated with non-coding SNVs, six exhibited significant upregulation, and five showed significant downregulation (Fig. 3d).

We next sought to determine whether these differentially expressed genes are involved in bile acid metabolism and transport, and thus contribute to the observed phenotypic changes. The bacterial bile acid detoxification pathway is now well-characterized as a coordinated process of stepwise modification and efflux. First, intracellular bile salt hydrolase (cbh) catalyzes the deconjugation of conjugated bile acids (e.g., glycocholic acid) into free bile acids and amino acids. The liberated bile acids can then be eliminated by one of two routes: (1) direct efflux via Major Facilitator Superfamily (MFS) efflux pumps (e.g., mdr, lmrB), which actively export bile acids across the cell membrane33; or (2) deeper modification followed by efflux, in which hydroxysteroid dehydrogenase (hdhA, 7β-HSDH) mediates oxidoreductive modification of free bile acids, and 7α-dehydroxylase enzymes encoded by the bai gene cluster subsequently catalyze their conversion into secondary bile acids. These secondary bile acids are subsequently extruded by the BaiF/BaiG transport system, thereby mitigating bile acid-induced membrane permeability and cytotoxicity34 (Fig. 3e). In total, fifteen key genes coordinate this detoxification cascade.

By cross-referencing the 59 (48 associated with coding SNVs and 11 associated with non-coding SNVs) differentially expressed genes with those involved in bile acid metabolism (n = 15), we identified two key candidates that were significantly upregulated in the adapted strain. The first, cbh (K01442, EC 3.5.1.24), encodes bile salt hydrolase, which catalyzes the hydrolysis of conjugated bile acids into free bile acids and amino acids (log₂FC = 0.4802, FDR < 0.001; Fig. 3f). This gene is associated with a non-coding SNV (rs1846941), suggesting a potential regulatory variant affecting its expression. The second gene mdr (via NCBI BLASTp) was also significantly upregulated (log₂FC = 0.5355, FDR < 0.001; Fig. 3g). Mdr encodes an MFS multidrug resistance (MDR) efflux pump responsible for exporting free bile acids. We found that a non-synonymous SNV (rs402499) within its coding sequence may alter transporter function. Together, integration of genomic and transcriptomic data from adapted and original strains potentially highlights these two genes, both genetically and transcriptionally remodeled, as drivers of enhanced bile acid metabolism during host-mediated in vivo induction strategy.

Validation of bile acid metabolic enhancement mediated by two key SNVs
To test whether the two key SNVs identified in the adapted strain indeed contribute to enhancing bile acid metabolism and associated phenotypic improvements, we performed systematic validation. Firstly, given the potential for false positives in WGS, we confirmed the presence of the two SNVs using an additional Sanger sequencing. Compared to the parental genome, we observed: (1) a C-to-A substitution at position 1,846,941 in the upstream non-coding region (likely the promoter) of cbh; and (2) an A-to-G substitution at position 402,499 in the coding region of mdr, resulting in a missense mutation that changed the 192nd amino acid of the mdr protein from asparagine (Asn) to aspartic acid (Asp, N192D). These results fully agree with our WGS findings (Fig. 4a, b; “Methods”).

We next sought to confirm that the upregulation of cbh in the adapted strain was attributable to the non-coding mutation. To this end, we constructed transcriptional reporters in which either the wild-type or the rs1846941-harboring upstream sequence was fused to an mCherry cassette and introduced these plasmids into the parental B. animalis subsp. lactis H22B656. Under three different bile-salt conditions (0.3% deoxycholic acid, DCA; cholic acid, CA; and chenodeoxycholic acid, CDCA), all strains carrying the rs1846941 variant exhibited significantly higher mCherry fluorescence than those with the wild-type sequence, indicating that this non-coding mutation enhances cbh transcription during bile acid stress (fold changes of 1.30, 1.71, and 12.36 for DCA, CA, and CDCA, respectively; Multiple Mann-Whitney tests with FDR correction, q < 0.01; Fig. 4c).

To directly verify its regulatory effect within the genomic context, we knocked out the cbh promoter region containing rs1846941 and measured cbh expression levels under three different bile salt conditions (0.3% DCA, CA, and CDCA). The results showed that under these three different bile salt conditions, the cbh expression in the wild-type adapted strain was significantly higher than in the parental wild-type strain, suggesting that this non-coding mutation enhances cbh transcription under bile acid stress (fold changes of 1.65, 1.51, and 1.49 for DCA, CA, and CDCA, respectively; one-way ANOVA with Tukey’s HSD test, p < 0.01; Fig. 4d). Importantly, after knockout of the promoter sequence, cbh expression was significantly reduced, further confirming that the upregulation of cbh in the adapted strain is indeed driven by this non-coding mutation. Moreover, assessment of bile salt hydrolysis using the ninhydrin colorimetric assay revealed that the adapted strain exhibited significantly higher hydrolytic activity toward glycine-conjugated bile salts compared to the parental strain (fold change = 1.25; one-way ANOVA with Tukey’s HSD test, p < 0.01; Fig. 4e). For taurine-conjugated bile salts, hydrolytic activity was also higher in the adapted strain (fold change = 1.16). Notably, knockout of the promoter region significantly reduced hydrolytic activity toward conjugated bile salts in both strains. Together, these findings suggest that the promoter region containing rs1846941 regulates the cbh gene expression and plays a critical role in enhancing the hydrolysis of conjugated bile acids.

For rs402499, we similarly performed a knockout of the mdr gene, which carries the mutation, and assessed mdr expression under three different bile salt conditions (0.3% DCA, CA, and CDCA) to determine whether the mutation directly affects gene expression. Under these three bile acid treatment conditions, the expression of mdr in the adapted strain was significantly higher compared to the parental wild-type strain, suggesting that this mutation enhances mdr transcription in response to bile acid stress (fold changes of 1.64, 2.13, and 1.73 for DCA, CA, and CDCA, respectively; one-way ANOVA with Tukey’s HSD post-hoc test, p < 0.01; Fig. 4f). Further analysis revealed that knockout of the mdr gene resulted in a significant reduction in mdr expression, confirming that the mutation directly impacts the phenotype through the regulation of mdr expression. Consistently, heterologous expression assays confirmed that the N192D mutation enhances mdr expression at translational levels under bile salt stress (Supplementary Fig. 2). Additionally, we performed complementation experiments to assess whether the restoration of mdr expression could recover the phenotype. Under CA, DCA, and CDCA treatment conditions, the complemented adapted strain exhibited significantly higher mdr expression than the complemented wild-type strain (One-way ANOVA with Tukey’s HSD post-hoc test, p < 0.01; Fig. 4f). These findings suggest that the rs402499 (N192D) mutation drives adaptive phenotypes by increasing mdr expression under bile acid stress.

MDR, as a transporter protein, facilitates the efflux of free bile acids, thereby mitigating bile acid toxicity and improving bacterial survival35. To further validate the causal relationship between mdr gene expression and the phenotype, we assessed the viable bacterial counts under bile acid stress conditions. The results showed that the colony counts of the adapted wild-type strain were significantly higher than those of the parental wild-type strain under all three bile salt conditions (One-way ANOVA with Tukey’s HSD post-hoc test, p < 0.01; Fig. 4g). Notably, under CA and CDCA treatment, the mdr knockout strains exhibited a marked growth disadvantage, whereas the differences between knockout strains were less pronounced under DCA conditions (p = 0.5499). In complementation experiments, the complemented adapted strain showed significantly higher colony counts than the complemented wild-type strain, further supporting the functional role of mdr gene expression in bile acid stress. Together, the N192D substitution significantly impacted both the transcriptional and translational levels of the mdr gene under bile stress, conferring a robust expression advantage and increasing the strain’s survival rate.

Since rs402499 encodes a non-synonymous N192D substitution, we hypothesized that this change might also alter the protein’s structural conformation in a way that improves its ability to export free bile acids. To test this hypothesis, we conducted comparative structural modeling and molecular docking simulations for both the wild-type and N192D variant of the protein to assess their binding affinity for free bile acids (“Methods”). The simulations revealed that the mutant protein exhibited significantly lower binding free energy with all three substrates (DCA, CA, and CDCA) and formed a greater number of hydrogen bonds (two-tailed Mann-Whitney test with FDR-corrected, q < 0.001; Fig. 4h, i). These results suggest that the N192D mutation stabilizes the protein-substrate complex, thereby potentially facilitating more efficient substrate recognition and transmembrane transport.

Adapted strain outperforms the original strain in NAFLD prevention
To assess the beneficial potential of B. animalis subsp. lactis W5S9 in modulating host metabolism, we employed a diet-induced nonalcoholic fatty liver disease (NAFLD) mouse model, grounded in the established link between bile acid dysregulation and NAFLD progression36. We compared the adapted strain (W5S9) with its original strain (H22B656) in safety, bile acid metabolism, and preventive efficacy (Fig. 5a). Specific pathogen-free (SPF) C57BL/6J mice (n = 6 per group) were randomized into four cohorts: (1) Control (normal diet, “Methods”); (2) HFAD (normal diet + 60% fat); (3) HFAD-H22B656 (HFAD with daily oral gavage of 10⁹ CFU B. animalis subsp. lactis H22B656); and (4) HFAD-W5S9 (HFAD with daily oral gavage of 10⁹ CFU B. animalis subsp. lactis W5S9). Probiotics were administered in 200 µL saline by gavage throughout eight weeks, and the control or HFAD group also received empty 200 µL saline to minimize the effect of stimulation by oral gavage. At the end of the trial, we collected serum for biochemical profiling, harvested livers for histopathology, and obtained fecal samples for metagenomic sequencing and targeted bile acid metabolomics. Metagenomic analysis confirmed that there was no cross-contamination in the control and HFAD groups, and that oral gavage efficiency was equivalent in the HFAD-H22B656 and HFAD-W5S9 groups, consistent with the fecal qPCR quantification results (Fig. 5b, c and Supplementary Fig. 3a).

We then compared the effects of the adapted and original strains on NAFLD progression across three parameters: body weight, physiological markers, and histopathologic changes. Over the eight-week intervention, all HFAD-fed groups (e.g., HFAD, HFAD-H22B656, and HFAD-W5S9) gained more weight than the control group on standard chow. Within the HFAD-fed groups, HFAD-W5S9 mice exhibited the smallest weight gain (two-tailed Mann-Whitney test, p = 0.0931), followed by HFAD-H22B656 (p = 0.0411), while the HFAD group showed the greatest increase (p = 0.0022). By week 8, body weight in the HFAD-W5S9 cohort was significantly lower than in the HFAD group (p = 0.0022), whereas there was no significant difference between HFAD-H22B656 and HFAD at week 8 (p = 0.1905; Fig. 5d). To quantify this therapeutic effect, we analyzed the net weight gain. W5S9 intervention significantly prevented HFAD-induced weight gain by 33.6% (mean gain: 7.76g vs. 11.68g in the HFAD group; Two-tailed Mann-Whitney test, p = 0.0022), exhibiting the smallest weight gain among the HFAD-fed groups. In contrast, the original strain showed a less pronounced reduction with higher variability (mean gain: 8.14g vs. 11.68g in the HFAD group; inhibition rate = 30.3%, p = 0.1320, Fig. 5e). No mortality occurred in mice receiving either H22B656 or W5S9 throughout the eight weeks. In vivo and in vitro safety and stability assessments further demonstrated that the adapted strain W5S9 retains parental-like safety and genetic stability without detectable off-target effects (Supplementary Figs. 4 and 5).

Regarding serum physiological metrics, both the original strain and the adapted strain significantly reduced levels of liver injury marker AST, low-density lipoprotein cholesterol (LDL-C), and pro-inflammatory cytokines (IL-6 and TNF-α) by week 8 compared to the HFAD group (Two-tailed Mann-Whitney test, HFAD vs. HFAD-H22B656 or HFAD vs. HFAD-W5S9, p < 0.05). Notably, compared to the HFAD-H22B656 group, the HFAD-W5S9 group exhibited further improvements in high-density lipoprotein cholesterol (HDL-C), malondialdehyde (MDA), and total cholesterol (TC). For example, while both strains modestly increased serum HDL-C levels, the increase was significantly greater in the HFAD-W5S9 group than in the HFAD-H22B656 group (p = 0.041). For MDA levels, H22B656 did not significantly reduce lipid peroxidation, whereas HFAD-W5S9 had a marked decrease (p = 0.016). A similar trend was observed for TC, with W5S9 showing a significant reduction compared to H22B656 (p = 0.041; Fig. 5f).

Histopathological analysis further corroborated the physiological improvements. Hematoxylin and eosin (H&E) staining of liver sections showed that HFAD-H22B656 mice had moderate alleviation of hepatic steatosis relative to HFAD alone, although scattered microvesicular lipid droplets persisted. In contrast, livers from the HFAD-W5S9 group largely maintained normal lobular architecture, exhibiting only mild, focal microvesicular steatosis (Fig. 5g). Oil Red O staining confirmed a marked reduction in hepatic lipid accumulation in HFAD-W5S9 mice, with more pronounced histological recovery than in the HFAD-H22B656 mice (Fig. 5h). Quantitative image analysis of lipid droplet area (n = 6 per group) further reinforced this observation, revealing a significantly greater reduction in hepatic lipid content in the HFAD-W5S9 group compared to the HFAD group (p = 0.0225, Kruskal-Wallis with Dunn’s test), approaching levels observed in the healthy control group (Fig. 5i). Together, these results demonstrate that the in vivo-induced and directed-evolved probiotic W5S9 provides significantly greater protection against NAFLD progression than the original strain.

Adapted strain’s impact on bile acid-metabolism and NAFLD prevention
Given the strong interplay between gut microbiota and bile acid metabolism34, we sought to determine how W5S9 enhances this interaction under HFAD conditions. We compared the effects of the original H22B656 and adapted W5S9 on gut microbial composition and intestinal bile acid metabolism using metagenomic sequencing and targeted bile acid profiling. In terms of taxonomic structure, significant increases in microbial α-diversity were observed in both HFAD-H22B656 and HFAD-W5S9 mice relative to HFAD alone, as reflected by higher richness and Shannon indices (two-tailed Mann–Whitney test, p < 0.05; Fig. 6a). For β-diversity, principal coordinate analysis (PCoA) of species-level Bray-Curtis dissimilarity showed that, after eight weeks, the gut microbiota of both HFAD-H22B656 and HFAD-W5S9 mice diverged markedly from those of the HFAD and control groups (Fig. 6b). However, there was substantial overlap between the microbiota profiles of the HFAD-H22B656 and HFAD-W5S9 mice, with no significant difference observed between them (p = 0.73; Fig. 6c). Moreover, among 423 species shared by the two groups, no significant differences in relative abundance were detected (multiple Mann–Whitney tests, q > 0.05; Supplementary Fig. 3e), indicating that the adapted strain exerts a broadly similar impact on community structure as the original strain, and that W5S9’s functional enhancement did not come at the expense of major shifts in gut microbiota (Supplementary Fig. 3b, c, d).

Given the minimal impact of W5S9 on overall taxonomic composition relative to its original strain, we next examined functional alterations, focusing on bile acid metabolism. We performed pairwise comparisons of targeted fecal bile acid profiles across the four groups using orthogonal partial least squares discriminant analysis (OPLS-DA). We found that: (1) HFAD and control groups showed clear separation, with 31 significantly altered bile acid metabolites (VIP > 1, |log₂FC| > 1; Fig. 6d); (2) HFAD and HFAD-H22B656 also diverged distinctly, identifying 28 differential metabolites (Fig. 6e); (3) HFAD-W5S9 profile deviated even further from HFAD, with 25 significantly altered bile acids. Notably, glycochenodeoxycholic acid (GCDCA) and taurochenodeoxycholic acid (TCDCA) were markedly reduced in B. animalis subsp. lactis W5S9-treated mice (Fig. 6f); (4) A clear separation between HFAD-H22B656 and HFAD-W5S9 emerged, with eight differential metabolites identified (Fig. 6g). Together, these metabolomic data demonstrate that W5S9 exerts a stronger modulatory effect on bile acid-related pathways than the original strain. Consistent with these metabolomic changes, HUMAnN3-based functional profiling revealed higher total bile salt hydrolase (EC 3.5.1.24) gene abundance in HFAD-W5S9 than HFAD-H22B656 group, with stratified analysis indicating that this enhancement is primarily driven by W5S9 colonization and enzymatic activity rather than broad remodeling of the resident microbiota (Supplementary Fig. 6).

We next mapped the significantly altered metabolites from our OPLS-DA comparisons onto the primary bile acid biosynthesis pathway (KEGG KO00120). To compare the metabolic modulation capabilities of the original and adapted strains, we generated a Venn diagram highlighting shared and unique pathway metabolites. Both B. animalis subsp. lactis H22B656 and W5S9 significantly altered levels of GCDCA, glycocholic acid (GCA), chenodeoxycholic acid (CDCA), and cholic acid (CA). Notably, B. animalis subsp. lactis W5S9 additionally regulated TCDCA and taurocholic acid (TCA), indicating broader engagement with primary bile acid biosynthesis (Fig. 6h). Furthermore, KEGG enrichment analysis revealed that differential metabolites in the HFAD vs. Control group comparison, as well as HFAD vs. HFAD-W5S9, were significantly enriched in KO00120 (p = 0.0357). In contrast, metabolites distinguishing HFAD from HFAD-H22B656 mice did not show significant enrichment (p = 0.2143; Fig. 6i).

In the primary bile acid biosynthesis pathway, primary bile acids are conjugated with glycine or taurine in the liver to yield TCDCA, GCDCA, TCA, and GCA, which are then secreted into the intestine via the biliary tract. Approximately 95% of these conjugated bile acids undergo active reabsorption by intestinal epithelial cells, while the remaining fraction is deconjugated by microbial BSHs, producing free bile acids such as CDCA and CA (Fig. 6j). To comprehensively evaluate the regulatory magnitude of the strains, we analyzed the log fold changes (log₂FC) relative to the HFAD model (Supplementary Fig. 7). While the parental strain H22B656 elicited a moderate reduction in pathological bile acid accumulation, the adapted strain W5S9 exhibited a significantly more profound restorative effect. Specifically, W5S9 treatment resulted in a substantial decrease in both conjugated bile acids (TCDCA, log₂FC = –1.14; TCA, log₂FC = –1.19; GCDCA, log₂FC = –∞; GCA, log₂FC = –2.68;) and free primary bile acids (e.g., CDCA, log₂FC = –2.09; CA, log₂FC = –2.02) in the intestinal lumen (Fig. 6k). Notably, bile acid concentrations in HFAD-W5S9 mice approached those of the Control group more closely than in the HFAD-H22B656 mice, indicating a more robust reestablishment of bile acid homeostasis. Given the central role of primary bile acid biosynthesis in maintaining enterohepatic equilibrium and regulating hepatic lipid metabolism, and its dysregulation as a key driver of NAFLD, these data suggest that the adapted strain enhances bile acid transformation and clearance, thereby helping to reestablish metabolic balance and mitigate disease progression.

Methods
Animal ethics and housing conditions
All animal experiments were approved by the Animal Ethics Committee of Hainan University (Approval No. HNUAUCC-2023-00190) and were conducted in strict accordance with the Guidelines for the Care and Use of Experimental Animals of Hainan University. Germ-free mice were maintained in sterile isolators, while specific pathogen-free (SPF) mice were housed in SPF animal facilities. All mice were kept under a 12 h light/12 h dark cycle and housed individually, with a controlled ambient temperature of 25 ± 2 °C and relative humidity of 55 ± 5%.

Calcium precipitation assay for bile salt hydrolase activity
Bifidobacterium animalis subsp. lactis (B. animalis subsp. lactis) strains were obtained from the Tropical Probiotics Development and Utilization Team at Hainan University, which were originally isolated from the intestinal tracts of healthy infants.

To screen for B. animalis subsp. lactis strains with bile acid-metabolizing capacity, 12 randomly selected strains were subjected to a calcium precipitation assay to evaluate their bile salt hydrolase activity. Briefly, 50 µL of each overnight culture was spotted onto MRS-cysteine (MRS-cys) agar supplemented with 0.1% cysteine·HCl, 3.0 g/L taurodeoxycholic acid sodium salt, 3.0 g/L glycodeoxycholic acid sodium salt, 2.0 g/L sodium thioglycolate, and 0.37 g/L anhydrous CaCl₂. Plates were incubated anaerobically at 37 °C for 48 h, after which the diameter of the calcium precipitation halo (indicative of bile salt hydrolase activity) was measured. MRS-cys agar without bile acids served as the negative control. All assays were performed in triplicate.

Bile salt hydrolases specific activity assay
To determine the BSH-specific activity of the candidate strains, we applied the classical ninhydrin colorimetric assay to quantify released amino acids55, a well-established method for amino acid measurement that indirectly reflects BSH hydrolysis efficiency. Specifically, overnight cultures of candidate strains were harvested at mid-log phase and thoroughly mixed. A 1 mL aliquot of culture was centrifuged at 3000 × g for 10 min at 4 °C, and the supernatant was discarded. The cell pellet was washed twice with 200 µL of ice-cold 0.1 M PBS (pH 6.0) and then resuspended in 500 µL of PBS (pH 6.0) to prepare a concentrated cell suspension.

BSH-specific activity was then assayed using glycine-conjugated bile acids (10 mM) as the substrate. In each reaction, 100 µL of the prepared cell suspension was mixed with 100 µL of substrate solution, and the final volume was adjusted to 1 mL with PBS (pH 6.0). Mixtures were incubated at 37 °C for 30 min. In parallel, blank controls contained 100 µL of cell suspension and 900 µL of PBS (no substrate). At the end of the incubation, 500 µL of 15% trichloroacetic acid was added to terminate the reaction, followed by 3 min of incubation at room temperature. Samples were then centrifuged at 14,000 × g for 10 min, and the supernatant was transferred to a 15 mL tube. To each supernatant, 1.5 mL of ninhydrin reagent was added and gently mixed before heating in a boiling water bath for 15 min. Tubes were immediately cooled on ice for 3 min, and absorbance was measured at 570 nm. Released amino acid concentrations were calculated against a concurrently prepared standard curve.

Simultaneously, total protein concentration in the crude lysate was determined using the Bradford assay with bovine serum albumin (BSA) as the standard. Specific activity (U/mg) was defined as the micromoles of amino acid released per milligram of protein per minute, and calculated using the following formula:

Specific activity (U·mg⁻¹) = Amount of amino acid produced (μmol) / (Protein content (mg) × Reaction time (min)) (1)

To assess functional enhancement in evolved strains relative to their parental counterparts, we calculated the fold change in BSH activity as follows:

Fold change = Cₙ / C₀ (2)

where Cₙ is the Specific activity by the nth induced strain, and C₀ is the corresponding value for its original strain (n = 1, 2, ..., 60).

Assessment of tolerance and growth performance
To comprehensively evaluate the gut-survival potential of these 12 candidate strains, we assessed their acid tolerance, bile acid tolerance, and growth performance. For acid and bile acid tolerance assays, strains were grown to mid-log phase and then transferred (2% v/v) into modified MRS-cys broth: one set adjusted to pH 2.0 (acid tolerance) and another supplemented with 0.3% (w/v) bovine bile (bile acid tolerance). Acid-treated cultures were incubated anaerobically at 37 °C for 2 h, and bile acid-treated cultures for 4 h, with untreated cultures in MRS-cys as controls. After incubation, OD₆₀₀ was measured in both treated and control groups, and survival rates (%) were calculated. As for growth performance, activated cultures were inoculated at 2% (v/v) into MRS-cys broth and incubated anaerobically at 37 °C for 36 h. OD₆₀₀ measurements were taken every 4 h, with uninoculated medium serving as the control.

In vitro induction experiment
To specifically enhance bile acid-metabolizing capacity, each candidate strain was first inoculated into 25 mL of MRS-cys broth and incubated anaerobically at 37 °C until reaching the logarithmic growth phase. At that point, 2% (v/v) of the culture was transferred into an equal volume of fresh MRS-cys medium supplemented with 0.3% (w/v) conjugated bile acids (a 1:1 mixture of sodium glycodeoxycholate and sodium taurodeoxycholate) and incubated anaerobically for 24 h. This serial passaging (using a 2% (v/v) inoculum each time) was repeated for ten consecutive rounds. After the tenth passage, the evolved population was diluted and plated onto solid MRS-cys agar. Subsequently, 60 colonies were randomly selected and assayed for BSH-specific activity to evaluate the effectiveness of the in vitro induction.

Host-mediated in vivo induction
To implement the in vivo adaptation strategy, we employed GF mice as an evolution bioreactor to eliminate interference from endogenous microbial strains. Male GF C57BL/6J mice (7 weeks old; n = 4) were obtained from Shenzhen Jingtai Biotech (Guangdong, China) and maintained in sterile isolators on irradiated standard chow for 7 days. To establish colonization by B. animalis subsp. lactis H22B656, each mouse received a 200 µL oral gavage of a 10⁹ CFU suspension once daily for three consecutive days. Following colonization, mice were subjected to a stepwise high-fat high-cholesterol dietary regimen over a 7-week induction period: weeks 1–3 on chow containing 30% fat and 0.2% cholesterol; weeks 3–5 on 45% fat and 0.2% cholesterol; and weeks 5–7 on 60% fat and 0.2% cholesterol.

To isolate adapted strains from GF mice, fresh stool samples were aseptically collected at the end of weeks 1, 3, 5, and 7, immediately suspended in 1.5 mL of sterile saline, and thoroughly homogenized under anaerobic conditions. Homogenates were serially diluted to 10⁻⁴ to 10⁻⁶ and plated on MRS-cys agar. After anaerobic incubation at 37 °C for 24 h, colonies from an appropriate dilution were randomly picked.

To ensure a representative evaluation of the adaptation process, 15 isolates were randomly selected from the plates of each sampling time point (Weeks 1, 3, 5, and 7), resulting in a total library of 60 induced isolates. These isolates were inoculated into MRS-cys liquid medium, labeled sequentially, and stored at –80 °C for subsequent analyses. Finally, 60 induced isolates were assayed alongside the original B. animalis subsp. lactis H22B656 for BSH-specific activity to evaluate the efficacy of the host-mediated in vivo induction.

Intervention diets
The normal diet (ND) was designed in accordance with established nutritional guidelines and comprised casein, corn starch, lard, vitamins, minerals, and fiber, providing a balanced mix of protein, fat, and carbohydrates. To generate an intestinal high-bile-acid environment in GF mice, we formulated a high-fat/high-cholesterol diet by increasing lard content to supply 30%, 45%, or 60% of energy from fat and supplementing with 0.2% cholesterol. For the NAFLD intervention study, a high-fat diet (HFAD) was similarly prepared from the normal diet to provide 60% energy from fat, without cholesterol supplementation. All diets were sterilized by γ irradiation and provided by Anhui Kuibu Shuyu Biotechnology Co., Ltd. Full dietary compositions are listed in Supplementary Data 2.

Genome assembly and annotation of the original strain
To obtain the complete reference genome of the original strain B. animalis subsp. lactis H22B656, genomic DNA was extracted from pure cultures of B. animalis subsp. lactis H22B656 using the QIAamp DNA Stool Mini Kit (Qiagen), following the manufacturer’s protocol. DNA concentration and integrity were assessed by NanoDrop spectrophotometry (A260/A280 ≈ 1.8) and 0.8% agarose gel electrophoresis. Samples meeting these quality criteria were sent to Beijing Novogene Bioinformatics Technology Co., Ltd. for third-generation long-read sequencing on the Oxford Nanopore Technologies (ONT) PromethION platform.

Nanopore libraries were prepared using the ONT Ligation Sequencing Kit (Oxford Nanopore Technologies) as follows: Genomic DNA integrity was assessed by gel electrophoresis and quantified using a Qubit fluorometer. The DNA underwent end-repair and A-tailing reactions, after which sequencing adapters and motor proteins were ligated to both ends. Unligated molecules were removed via purification steps, and libraries meeting quality standards were sequenced on the PromethION platform. Sample B. animalis subsp. lactis H22B656 generated 433,564 raw reads with a total output of 2.4 Gb.

For genome assembly, raw reads were first demultiplexed and adapter-trimmed to generate subreads, which were then filtered to retain only those ≥ 50 bp in length. De novo assembly of the filtered subreads was performed with Canu56 (v2.2). The resulting draft contigs were polished through three rounds of alignment and consensus correction by coupling Minimap257 (v2.22) and Racon58 (v1.4.13, command: minimap2 -t 40 assembly.fasta pass.fastq.gz > round1.paf; racon -t 4 pass.fastq.gz round1.paf assembly.fasta > round1.fasta).

Gene prediction and functional annotation of the complete genome were performed with Prokka59 v1.14.6 (command: prokka --outdir test --prefix test test.fa). Predicted proteins were further annotated via the KEGG Automatic Annotation Server60 (KAAS v4.0; http://www.genome.jp/tools/kaas/), using the bi-directional best hit method for pathway assignment. To further annotate hypothetical proteins, BLASTp searches were also conducted with NCBI BLAST+ v2.16.0 against the non-redundant (nr) protein database (e-value ≤ 1e-5, identity ≥ 30%, query coverage ≥ 70%, up to 20 alignments; command: blastp -query query.faa -db nr -out blast.tsv -outfmt 6 -evalue 1e-5 -num_threads 24 -num_alignments 20).

Whole genome sequencing of the adapted strains
To investigate the genetic basis of its enhanced bile acid-metabolizing capacity, we performed whole-genome sequencing (WGS) of the adapted strain W5S9. Genomic DNA from purified W5S9 cultures was extracted and assessed for its purity (A260/A280 = 1.8-2.0) and fragment integrity (>20 kb). DNA was then randomly sheared to the desired fragment size using a Covaris ultrasonicator. Fragmented DNA underwent end repair and 3’ A-tailing, followed by ligation of Illumina sequencing adapters. Adapter-ligated libraries were size-selected to remove off-target fragments, PCR-amplified to enrich adapter-containing molecules, and purified to eliminate residual contaminants. Equimolar libraries were pooled and sequenced by Novogene Co., Ltd. (Beijing, China) using an Illumina HiSeq 2500 platform (2 × 150 bp paired-end). Raw reads were stringently filtered for quality using fastp61 v0.23.2. Sample W5S9 yielded 8,310,066 raw reads, of which 8,289,624 high-quality reads were retained after filtering.

SNVs calling
To accurately capture genetic variants acquired during in vivo adaptation, we performed SNV calling on WGS data using two complementary pipelines (inStrain30 v1.0.0 and BCFtools31 v 1.21) and retained only the intersection of their call sets for downstream analysis. First, WGS data were aligned to the original H22B656 reference genome with Bowtie262 v2.5.4. Alignments were converted to coordinate-sorted, indexed BAM files using SAMtools31 v1.22. SNV detection was then carried out independently with: (1) inStrain (command: inStrain profile sample.bam ref.fna -c 100 -f 0.49 -o sample.profile -p 24 -g ref_genes.fna); (2) BCFtools (command: bcftools mpileup -C 50 -Ou -m 3 -F 0.0002 -f ref.fa sorted.bam | bcftools call -mv -Ob -o variants.bcf). High-confidence SNVs were defined as variants identified by both methods. We identified a total of 282 high-confidence SNVs in the adapted strain W5S9. The distribution of SNVs across the genome was plotted using the R package circlize (version 0.4.15, Fig. 3a).

Promoter prediction for non-coding SNVs
To determine whether SNVs in non-coding regions might influence downstream gene expression, we conducted promoter prediction and target-gene association analyses for all intergenic SNVs. First, each SNV’s nearest gene was identified using BEDTools63 v2.18, and the distance to the gene’s transcription start site (TSS) was calculated. Genes located on the same strand, with SNVs positioned within 0–300 bp upstream of their TSS, were selected as candidate regulatory targets. Promoter regions encompassing each candidate SNV were then predicted using BPROM (www.softberry.com), and the predicted promoter elements were inspected for canonical −10/−35 motifs to evaluate their regulatory potential.

Transcriptome sequencing and data analysis
To identify key genes with expression changes caused by genetic variants acquired during in vivo adaptation, we performed transcriptome sequencing (RNA-seq) on the original strain H22B656 and the adapted strain W5S9. Briefly, each strain was cultured anaerobically for 12 h in MRS-cys broth supplemented with 0.3% conjugated bile acid (sodium glycodeoxycholate and sodium taurodeoxycholate). Cultures were centrifuged at 1500 × g for 30 min, and pellets were immediately processed for RNA extraction.

Total RNA was isolated using a Qiagen RNeasy kit, and integrity (RIN > 8) and purity (A260/A280: 1.8-2.0) were confirmed by agarose gel electrophoresis and NanoDrop spectrophotometry. Ribosomal RNA was depleted, and the remaining RNA was purified by ethanol precipitation. First-strand cDNA synthesis was primed with random hexamers, and second-strand synthesis incorporated dUTP to preserve strand specificity. The resulting double-stranded cDNA underwent end repair, 3’ A-tailing, adapter ligation, size selection, USER enzyme digestion to remove the dUTP-containing strand, PCR enrichment, and bead-based purification to yield strand-specific libraries. Library concentration was measured by Qubit and qPCR, and fragment size distribution was assessed on an Agilent Bioanalyzer. Libraries meeting quality thresholds were pooled equimolarly and sequenced (2 × 150 bp) on an Illumina HiSeq 2500 by Novogene Co., Ltd. (Beijing, China).

Raw data were quality-filtered with fastp v0.23.2 to remove adapter contamination, trim low-quality bases (Q < 20), and discard reads containing ambiguous nucleotides (N), resulting in high-quality clean reads for downstream analysis. Sequencing of six transcriptomics samples yielded an average of 11.95 ± 0.15 million raw reads and 11.79 ± 0.15 million clean reads per sample (mean ± SEM). Transcript-level quantification was performed with Salmon64 v1.9.0 against the reference transcriptome (salmon quant -i mg -l A -1 test_1.fastq.gz -2 test_2.fastq.gz -o test.quant). Differential expression analysis between the original and adapted strains was performed in R (v4.4.1) using the DESeq2 (v1.44.0) package (n = 3 per strain, Benjamini–Hochberg FDR correction; FDR < 0.001 and |log₂FC| > 0.4), and the resulting differential features were visualized as volcano plots using OmicShare (https://www.omicshare.com/; Fig. 3c, d). Gene expression differences between adapted and original strains for cbh and mdr were shown as bar charts in GraphPad Prism v10 (Fig. 3f, g).

Validation of identified SNVs by targeted Sanger sequencing
To validate the candidate SNVs (rs1846941 and rs402499) identified by WGS, we designed locus-specific primer pairs (mrd-F/mrd-R and rs1846941-f/ rs1846941-r; Supplementary Data 3) to amplify 200 bp regions centered on each variant (100 bp upstream and downstream). Purified PCR products were submitted to Beijing Qinke Biotechnology Co., Ltd. for Sanger sequencing. Following sequencing, the raw reads were aligned against the B. animalis subsp. lactis H22B656 reference genome to verify the target mutations. Sequence visualization and variant confirmation were performed using Geneious Prime v2022.2 (Fig. 4a, b).

Determination of the genetic stability of key SNVs
To verify the genetic stability of key SNVs in the adapted strain, we conducted a rigorous long-term continuous passage experiment followed by whole genome sequencing (WGS). Specifically, the adapted strain W5S9 was subjected to daily serial passaging for 4 weeks under two conditions: (i) a non-selective medium (MRS broth supplemented with L-cysteine hydrochloride), used to assess stability in the absence of selection pressure; and (ii) a selective medium (MRS broth supplemented with bile salts and L-cysteine hydrochloride), used to mimic gut environmental pressure. Target strains were re-isolated weekly, genomic DNA was extracted, and WGS was performed to monitor potential genetic changes (Supplementary Fig. 4).

Reporter plasmid construction and fluorescence assays
To assess the regulatory impact of SNVs rs1846941 and rs402499 on target gene expression, we constructed mCherry and EGFP reporter plasmids harboring either the wild-type or variant sequences.

For SNV rs1846941 located in the non-coding region, we PCR-amplified 100 bp fragments (50 bp upstream and downstream) carrying either the wild-type or rs1846941 allele from genomic DNA of the original H22B656 and the adapted W5S9 strains, respectively. Each non-coding fragment was fused to mCherry by overlap-extension PCR, and the resulting expression cassettes were cloned into the EcoRI site of pLH01. Recombinant plasmids were introduced into B. animalis subsp. lactis H22B656 by electroporation. Single colonies were then grown anaerobically at 37 °C in MRS-cys medium (supplemented with 2% glucose, 10 µg/mL chloramphenicol, and 0.3% (w/v) conjugated sodium deoxycholate, sodium cholate, or sodium chenodeoxycholate) to an OD₆₀₀ of ~0.5. Cells were then harvested, washed twice with PBS, and fluorescence was measured on a BioTek Synergy H1 multimode plate reader.

For SNV rs402499 located in the coding region (in mdr), wild-type and N192D gene-egfp fusions were assembled by overlap-extension PCR and cloned into SalI/XhoI–linearized pET28a. Plasmids were transformed into E. coli BL21 (DE3) by heat shock. Transformants were cultured in LB medium (with 0.3% conjugated sodium deoxycholate, conjugated sodium cholate, or conjugated sodium chenodeoxycholate and 50 µg/mL kanamycin) at 37 °C to OD₆₀₀ = 0.6–0.8. After harvesting and two PBS washes, fluorescence was recorded on the Synergy H1. Statistical analyses and data visualization were performed using GraphPad Prism v10. For Fig. 4c and Supplementary Fig. 2, multiple two-tailed Mann–Whitney tests were used to compare H22B656 and W5S9 (n = 3 independent samples per group), and multiple comparisons were controlled using a two-stage step-up false discovery rate (FDR) procedure with Q = 1%.

Construction of gene knockout and complemented strains
To establish a causal relationship between the rs1846941 and rs402499 mutations and the phenotype of interest, we employed a homologous recombination strategy using the temperature-sensitive suicide plasmid pKO403 to generate deletion mutants targeting the mdr and cbh promoter regions in both the parental strain Bifidobacterium animalis subsp. lactis H22B656 and the adapted strain B. animalis subsp. lactis W5S9. Genomic DNA from the respective strains was used as the template to amplify the upstream and downstream homologous arms flanking the mdr gene and cbh promoter regions by PCR using primers listed in Supplementary Data 3. Amplicons were resolved by agarose gel electrophoresis and purified. The pKO403 plasmid was linearized by double digestion with SalI and BamHI and seamlessly assembled with the purified homologous arms using 2× Ezmax® Universal CloneMix (37 °C, 30 min). Assembly products were transformed into E. coli DH5α competent cells, and recombinant plasmids were identified by colony PCR. Plasmids confirmed by PCR were isolated and electroporated into electrocompetent H22B656 and W5S9 cells (2.1 kV, 200 Ω, 25 µF). Following electroporation, cells were recovered in 1 mL MRS-cys medium under anaerobic conditions at 37 °C for 4 h and plated on MRS-cys agar containing chloramphenicol (10 µg/mL) to select for single-crossover integrants. Chloramphenicol-resistant colonies were cultured in MRS-cys broth supplemented with chloramphenicol (10 µg/mL) at 37 °C for 48 h and subsequently plated on MRS-cys agar with chloramphenicol for anaerobic incubation at 42 °C for 48 h to promote double-crossover recombination. Successful deletion of the target region was confirmed by PCR, and double-crossover deletion mutants were obtained. To eliminate the pKO403 plasmid, verified mutants were passaged for three consecutive generations in antibiotic-free MRS-cys medium.

For complementation of the Δmdr mutant, a 1648-bp fragment containing the full-length mdr gene and its native promoter was PCR-amplified from the genomic DNA of H22B656 or W5S9 using primers listed in Supplementary Data 3 and cloned into SalI/BamHI-digested pKO403. The resulting plasmid was introduced into the Δmdr strain by electroporation, and transformants were selected on MRS-cys agar containing chloramphenicol (10 µg/mL). Successful construction of complemented strains was confirmed by PCR.

RT-qPCR gene expression analysis
To assess the impact of the rs402499 and rs1846941 variants on gene expression, bacterial cultures were grown in MRS-cys medium supplemented with 0.3% (w/v) conjugated sodium deoxycholate, sodium cholate, or sodium chenodeoxycholate. After 12 h of incubation, bacterial cells were collected for subsequent analysis. Total RNA was extracted with TRIzol reagent and treated with DNase I to remove genomic DNA. cDNA synthesis was performed using the PrimeScript™ RT kit according to the manufacturer’s instructions. The cycle threshold (Cₜ) values were normalized to the reference gene rplx to obtain ΔCₜ values. Relative gene expression was calculated using the comparative Cₜ method (2⁻ΔΔCT), with the wild-type Bifidobacterium animalis subsp. lactis H22B656 strain serving as the calibrator (relative expression = 1) for each bile salt condition independently. Statistical analyses were performed using R v4.4.1 software with the rstatix (v0.7.2) and ggpubr (v0.6.0) packages. Data are presented as mean ± SD, with individual data points from biological replicates overlaid (n = 3). Comparisons between strains under each bile salt condition were analyzed using one-way analysis of variance (ANOVA) followed by Tukey’s Honest Significant Difference (HSD) post-hoc test for multiple comparisons. Statistical significance was defined as p < 0.05 (Fig. 4d, f).

Bile acid stress survival assay for mdr phenotypic validation
To validate the contribution of mdr expression to bile acid tolerance, survival under bile acid stress was assessed by viable colony counting. The parental wild-type strain H22B656, the adapted wild-type strain W5S9, their isogenic mdr knockout mutants, and the corresponding complemented strains were grown overnight under standard culture conditions as described above. Cultures were harvested and adjusted to the same optical density, and equal aliquots were inoculated into fresh medium supplemented with 0.3% (w/v) cholic acid (CA), deoxycholic acid (DCA), or chenodeoxycholic acid (CDCA). After incubation under bile acid stress, cultures were serially diluted and plated on bile acid–free agar. Colony-forming units (CFU) were counted after incubation, and survival was expressed as CFU per mL. Each condition was tested in three independent biological replicates, and data were analyzed by one-way ANOVA followed by Tukey’s HSD post hoc test (Fig. 4g).

Molecular docking analysis
To investigate the impact of the nonsynonymous variant rs402499 on MDR structure and its binding affinity for free bile acids, we used the AlphaFold server (https://alphafoldserver.com) to predict the three-dimensional structures of both wild-type and mutant MDR from their amino acid sequences, exporting the resulting models as PDB files for downstream analyses. The three-dimensional structures of deoxycholic acid (CID: 222528, https://pubchem.ncbi.nlm.nih.gov/compound/Deoxycholic-acid), cholic acid (CID: 221493, https://pubchem.ncbi.nlm.nih.gov/compound/Cholic-acid), and chenodeoxycholic acid (CID: 10133, https://pubchem.ncbi.nlm.nih.gov/compound/Chenodeoxycholic-acid) were retrieved from PubChem, converted from SDF to PDBQT format with Open Babel (v3.1.1), and assigned Gasteiger charges.

Both wild-type and mutant MDR models were prepared for docking with AutoDockTools (v1.5.7) by removing crystallographic water molecules, adding hydrogens, and assigning Kollman charges before exporting them as PDBQT files. Potential ligand-binding pockets were identified using DoGSiteScorer (https://proteins.plus/). A ligand-centered docking grid was then defined (box size: 126³ Å³; grid spacing: 0.368 Å) and centered at coordinates x = -0.855, y = 10.657, z = -4.238 (dimensions: 20.33 × 20.33 × 20.33 Å³). Targeted molecular docking was performed with AutoDock Vina65 (v1.5.7). Statistical analyses and data visualization were performed using GraphPad Prism v10 (Fig. 4h). The binding pose with the lowest predicted binding energy was visualized in PyMOL (v2.5.4), and key ligand-receptor interactions, including hydrogen bonding, hydrophobic contacts, and shape complementarity, were analyzed (Fig. 4i).

Diet-induced NAFLD mouse model and probiotic gavage
To assess the beneficial potential of the adapted B. animalis subsp. lactis W5S9 in modulating host metabolism, we employed a diet-induced NAFLD mouse model in male SPF C57BL/6J mice (7 weeks old), based on the established link between bile acid dysregulation and NAFLD progression66. Mice were randomized into four groups (n = 6 per group): (1) Control (n = 6, normal diet, daily gavage with 200 µL sterile saline); (2) HFAD (n = 6, high-fat diet, daily gavage with 200 µL sterile saline); (3) HFAD-H22B656 (n = 6, high-fat diet, daily gavage with 200 µL sterile saline containing 10⁹ CFU of original strain B. animalis subsp. lactis H22B656); and (4) HFAD-W5S9 (n = 6, high-fat diet, daily gavage with 200 µL sterile saline containing 10⁹ CFU of adapted strain B. animalis subsp. lactis W5S9). The experiment lasted for 8 weeks. Body weight was measured weekly. At the end of the study, serum and liver tissues were collected for biochemical measurements and histopathological analysis to assess the preventive effects of the adapted and original strains on NAFLD. Additionally, fecal samples were collected for metagenomic sequencing and targeted bile acid metabolomics analysis.

Physiological and histopathological assessments
For the physiological markers assay, the levels of high-density lipoprotein (HDL-C), low-density lipoprotein (LDL-C), total triglycerides (TG), total cholesterol (TC), aspartate aminotransferase (AST), malondialdehyde (MDA), IL-6 and TNF-α in serum were measured according to the instructions of the ELISA kit (Shanghai Xinyu Biotechnology Co., Ltd.).

For histopathological analysis, liver specimens were rinsed in ice-cold sterile saline, and portions were either fixed in 4% paraformaldehyde or snap-frozen in liquid nitrogen. Sample processing, including paraffin embedding, sectioning, staining, and image capture, was performed by Wuhan Sewell Biotechnology Co., Ltd. (China). Semi-quantitative analysis of lipid droplet content was performed on Oil Red O-stained sections using Fiji (v2.17.0). Data visualization and statistical analyses were conducted in GraphPad Prism v10 using a Kruskal–Wallis test followed by Dunn’s multiple comparisons test (n = 6 mice per group, Fig. 5i).

Whole metagenome sequencing and data analysis
To compare the effects of the adapted versus original strains on the gut microbiota composition in the diet-induced NAFLD mouse model, fresh stool samples were collected at the end of week 8 and immediately stored at −80 °C for whole metagenome sequencing (WMS). Genomic DNA was extracted using the QIAamp DNA Stool Mini Kit (Qiagen) and evaluated for integrity by 0.8% agarose gel electrophoresis, purity via A260/A280 ratios on a NanoDrop spectrophotometer, and concentration with a Qubit 2.0 fluorometer. Paired-end libraries (2 × 150 bp) were prepared and sequenced on an Illumina HiSeq 2500 platform at Beijing Novogene Bioinformatics Technology Co., Ltd. Raw reads were quality-filtered using fastp v0.23.2 to remove adapter sequences and exclude read pairs with more than 50% of bases at Q ≤ 5 or over 10% ambiguous nucleotides (N bases). On average, each sample yielded 45.28 ± 0.56 million raw reads and 44.08 ± 0.57 million clean reads (mean ± SEM). Taxonomic profiling of the WMS data was performed using MetaPhlAn (v4.0.6)67 to identify microbial species and estimate their relative abundances. To specifically distinguish the direct enzymatic contribution of the supplemented strains from potential regulatory effects on the resident microbiota, we performed functional profiling of the metagenomic data using HUMAnN (v3.9) with the uniref90_koregrouping option to quantify the abundance of bile salt hydrolase (EC 3.5.1.24)–encoding genes (Supplementary Fig. 6).

qPCR quantification of B. animalis in fecal samples
To quantitatively assess the abundance of the original strain H22B656 and the adapted strain W5S9 in fecal samples, we designed H22B656-specific primers using the NCBI Primer-BLAST tool (H22BF: ACGGCACGATTCATGTGTCT; H22BR: CCAGATTCCATCCTCCTCGC). Fecal samples (0.1 g) were then subjected to DNA extraction using the TIANGEN Bacterial Genomic DNA Extraction Kit, and quantitative real-time PCR (qPCR) was performed with the above primers. qPCR data were analyzed and visualized using GraphPad Prism v10, and statistical significance was evaluated using the Kruskal–Wallis test followed by Dunn’s multiple comparisons test (Fig. 5c).

Targeted metabolomic profiling and data analysis of fecal bile acids
A total of 24 fecal samples (n = 6 biological replicates per group) were analyzed. Fecal bile acid extraction was performed using HPLC-grade acetonitrile and methanol (Merck) as solvents. Analytical standards (CNW, IsoReag, and Sigma-Aldrich) were prepared as 1 mg/mL stock solutions in methanol and stored at −20 °C. Working solutions were diluted to appropriate concentrations for calibration curves prior to analysis. Approximately 20 mg of ball-milled fecal sample was homogenized in 495 µL of methanol, spiked with 5 µL of internal standard solution (10 µg/mL), and subjected to protein precipitation at −20 °C. Following centrifugation at 13,200 × g for 10 min (4 °C) and filtration through a protein precipitation plate, supernatants were analyzed. To monitor instrument stability, quality control (QC) samples were prepared by pooling equal aliquots of all supernatants and analyzed every 10 samples.

Chromatographic separation was performed on an ExionLC™ AD UHPLC system (AB Sciex) equipped with a Waters HSS T3 C18 column (100 × 2.1 mm, 1.8 µm) maintained at 40 °C. Mobile phases consisted of (A) 0.01% acetic acid with 5 mmol/L ammonium acetate in water and (B) 0.01% acetic acid in acetonitrile. A linear gradient elution was performed at a flow rate of 0.35 mL/min: 5% to 40% B in 0–1 min; 40% to 50% B in 1–7 min; 50% to 75% B in 7–12 min; 75% to 95% B in 12–14 min; maintained at 95% B for 2 min (14–16 min); and re-equilibrated at 5% B from 16 to 17.5 min.

Mass spectrometric detection was conducted on an AB Sciex 6500+ QTRAP system coupled to an ExionLC™ AD UHPLC via electrospray ionization (ESI) in negative mode. Key parameters included: ion source temperature 550 °C, spray voltage −4500 V, and curtain gas 35 psi. Data acquisition employed scheduled multiple reaction monitoring (MRM) with Analyst software (v1.6.3), and quantification was performed using MultiQuant software (v3.0.3). Metabolites were identified by matching retention times and MS/MS spectra with the MetWare Database (MWDB) and authentic standards. Quantification was performed using MultiQuant software (v3.0.3) based on linear calibration curves. Fecal bile acid profiling was performed by MetWare (http://www.metware.cn/) based on the AB Sciex QTRAP 6500 LC-MS/MS platform. A detailed list of MRM transitions, collision energies, retention times, and mass spectra for all targeted bile acids is provided in Supplementary Data 4.

For targeted fecal bile acid metabolomics, orthogonal partial least squares discriminant analysis (OPLS-DA) was performed using the R (v4.4.1) package ropls (v1.36.0) for pairwise comparisons between the indicated groups (n = 6 per group). Differential metabolites were defined as those with a variable importance in projection (VIP) > 1 and |log₂ fold change| > 1 between the indicated groups. Score plots and heatmaps were generated with ggplot2 (v3.5.2) and pheatmap (v1.0.12), respectively (Fig. 6d–g). Intersection analysis of differential metabolites was conducted with the Evenn online tool (https://www.bic.ac.cn/test/venn/#/) and visualized as a Venn diagram (Fig. 6h). KEGG pathway enrichment results were displayed as bubble plots using ggplot2 (v3.5.2) (Fig. 6i), and quantitative data for key differential metabolites were plotted as bar graphs in GraphPad Prism v10 (Fig. 6k).

Statistics and reproducibility
Statistical analyses were performed using GraphPad Prism v10 and R software (v4.4.1). Data are presented as mean ± standard deviation (SD) or mean ± standard error of the mean (SEM), as indicated in the figure legends. Data distribution was assessed for normality. Comparisons between two groups were conducted using unpaired two-tailed Student’s t-tests or Mann–Whitney U-tests. For comparisons among more than two groups, one-way analysis of variance (ANOVA) followed by Tukey’s post-hoc test or Kruskal–Wallis tests followed by Dunn’s multiple comparisons test were applied. To correct for multiple testing in large-scale datasets (e.g., metabolomics and transcriptomics), the False Discovery Rate (FDR) was controlled using the Benjamini–Krieger–Yekutieli two-stage step-up method (or the Benjamini–Hochberg method for RNA-seq).

Specific R packages were employed for data visualization and analysis, including ggplot2 (v3.5.2) for general plotting, ggradar (v0.2) for radar plots, pheatmap (v1.0.12) for heatmaps, ggsignif (v0.6.4) for significance annotation, and vegan (v2.6-10) for diversity analysis. All schematic diagrams were created using BioRender (Zhang, J. (2026) https://BioRender.com/3qyx3m3). All experiments were performed with at least three independent biological replicates to ensure reproducibility.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this study, they propose an in vivo, host-mediated adaptation strategy that emulates natural intestinal selective pressures to direct probiotic evolution. Specifically, by exposing B. animalis subsp. lactis to a dynamically shifting bile acid milieu, which was induced through a high-fat, high-cholesterol diet, they isolated an adapted strain with markedly enhanced bile acid-metabolizing capacity. WGS and transcriptome analyses identified two key SNVs underpinning this phenotype: one variant elevated deconjugation of bile acids, while the other augmented efflux of free bile acids. Subsequently, in a diet-induced NAFLD mouse model, dietary supplementation with the adapted strain led to superior bile acid turnover, reduced inflammatory markers, and significantly greater protection against disease progression compared to the parental strain. Congratulations on your valuable study. After correcting all of the following, your manuscript can be accepted.\n1. A flowchart is needed to illustrate the process of selecting H22B656.\n2.The feces of model animals should be quantitatively tested (RT-qPCR? or other methods...) to reflect the stability of bacterial colonization.\n3.The X-axis of Figure 3c? The Log2(FC) ?\n4.Is there any data supporting the impact of H22B656 on the intestinal wall of animals (The level of IL-6?10?TNF-a?)?\n5.At the final time point, were intestinal segments tissue responsible for metabolizing secondary bile acids collected？And to assess the activation status of relevant targets and the network between bacteria, fermented substances, and the host?\n6.Line 750: The Metaphlan4?\n7.Figure6: The authors only administered H22B656 to germ-free mice. Why is the Shannon diversity index so high? Please provide specific changes in the intestinal bacterial composition.\n8.The heatmap of 7d-g were scale of row or col?\n9.The specific software parameters should be provided.", "reviewer_2": "Recommendation: Major Revision\nThis manuscript presents a novel and compelling approach to probiotic optimization through host-mediated directed evolution. The concept of using the gut as a living bioreactor to select for functionally enhanced probiotics is innovative and holds significant translational potential, particularly in the context of NAFLD. The multi-omics integration and functional validation of key mutations are strengths. However, several major concerns must be addressed before the manuscript can be considered for publication.\nMajor Concerns:\n1.Lack of mechanistic causality between mutations and phenotype. While the two SNVs (rs1846941 and rs402499) are associated with enhanced bile acid metabolism, direct evidence of causality is limited. Perform complementation or knockout experiments in the parental strain to confirm that these mutations alone are sufficient to confer the enhanced phenotype.\n2.Insufficient characterization of the adapted strain’s ecological impact. The metagenomic data show no major shifts in community structure between HFAD-H22B656 and HFAD-WSS9 groups, but functional changes in the microbiome (e.g., BSH activity of other taxa) are not assessed. Include functional metagenomics or metatranscriptomics to evaluate whether WSS9 influences the bile acid metabolism of other gut microbes.\n3.Limited assessment of long-term stability and safety. The adapted strain was evaluated in an 8-week intervention, but the genetic stability of the mutations under prolonged colonization or in complex microbiota is not tested. Include a stability assay (e.g., serial passaging in vitro or in vivo without selection pressure) to assess whether the mutations revert or confer unintended traits.\n4.Incomplete translation to human relevance. The study relies solely on a mouse model. The applicability of this approach to human microbiomes remains speculative. Discuss limitations and propose future steps, such as using humanized gnotobiotic models or organoids, to enhance translational relevance.\n5.Statistical and methodological clarifications needed. Some figure captions and methods lack detail (e.g., n-values, statistical tests used in metabolomics, multiple testing correction in RNA-seq). Provide full statistical details in figure legends and methods. Ensure that all comparisons are appropriately corrected for multiple testing.\nMinor Concerns:\n1.Figure 3e: The schematic of bile acid metabolism is helpful but could be better annotated with gene names and pathways.\n2.The discussion section could be more concise and focused on the broader implications of the work beyond NAFLD.\n3.Please clarify whether the adapted strain was tested for any off-target effects (e.g., growth kinetics, biofilm production, antibiotic resistance, virulence factors, or other traits).", "reviewer_3": "This is a very interesting topic with the potential for application to address issues in human health and disease. The main premise of the study is to determine whether the “host gut … could serve as a bioreactor for the directed evolution of probiotic strains”. The authors utilize a probiotic strain of Bifidobacterium animalis with intrinsic bile acid-metabolizing capacity. When an in vivo approach was taken in germ-free mice under bile-acid stress, the probiotic strain exhibited an increase in bile acid metabolism whereas in vitro evolution methods failed to produce any functional gains. The authors identified 2 single nucleotide variants that could explain the enhanced bile acid metabolism in the adapted strain. They go on to show the functional relevance of the adapted strains in a diet-induced mouse model of nonalcoholic fatty liver disease (NAFLD).\nThe noteworthy results include the efficacy of the in vivo induction strategy for enhancing bile acid-metabolism, the identification and functional validation of two SNVs that potentially confer the enhanced bile acid metabolism. The data presented with respect to the two SNVs is for the most part convincing. However, effects of the two SNVs on enhanced bile acid metabolism seems over-interpreted in some sections and needs clarification in other areas. I think where the findings could also be strengthened is in the studies focusing on the translational potential in the NAFLD model.\nMain comments:\nFigure 4d: The authors state that “mdr expression was significantly higher in the mdr-N192D under all bile acid conditions”. While that might be accurate at the 3h time point, it does not fully describe the data shown about mdr expression going down from baseline for both WT and N192D variant following induction with cholic acid CA and chenodeoxycholic acid CDCA. In addition, the followup sentence “Together, these results demonstrate that the N192D substitution upregulates mdr expression” isn’t entirely accurate…could this statement be revised given that the variant upregulates mdr expression to a greater degree after induction specifically with deoxycholic acid induction and not with CA or CDCA? The second half of the statement “potentially by altering local mRNA secondary structure (thus enhancing ribosome binding or translation initiation), optimizing codon usage bias, modulating mRNA folding dynamics or modifying regulatory elements affecting splicing translation” seems overly speculative.\nLine items 272-273: in the results section, it states “N192D mutation stabilizes the protein-substrate complex, thereby facilitating more efficient substrate recognition and transmembrane transport” is declarative without evidence…would it be more accurate to state “…complex, potentially facilitating more efficient substrate recognition…”?\nAdapted strain on weight gain: The three parameters chosen are useful outcome measurements to gauge the effect of the adapted strains. However, what could be clarified is the main comparison of interest. In the weight gain figure (Figure 5c), the authors state the HFAD-W5S9 treated mice showed the smallest weight gain from control, followed by the parental strain treated mice. The effect size of the adapted strain is not immediately accessible to the reader. What is the absolute difference in effect?…did it prevent weight gain imposed by a HFD by 30%, 50%, etc. and how does that compare to the parental strain? The results section could be rewritten to clarify this point.\nAdapted strain on liver histology: It would be helpful to perform semi-quantitative scoring for histology and oil red O staining across samples within each group to get a more accurate assessment of the adapted strain’s impact, if any, on the liver. This is critical as non-alcoholic steatohepatitis rests ultimately on histologic diagnosis. Both the parent and adapted strain appeared to reduce AST levels so it is unclear if there are significant differences in their impact on liver pathology.\nFigure 6 figure legend title “The adapted B. animalis subsp. Lactis W5S9 re-balances intestinal bile-acid metabolism and prevents NAFLD progression” needs to re-worded, as “prevents NAFLD progression” is not shown in this figure.\nAdapted strain on bile acid metabolites: In figure 6K, the authors note that bile acid concentrations in the adapted strain treated mice approached those of the control group more closely than in the parental strain treated mice, “indicating a more robust reestablishment of bile acid homeostasis”. This reviewer finds it interesting that the probiotic parental strain also impacted bile acid metabolism. The data is presented in such a way that it is difficult to gauge the effect of parental strain vs the adapted strain in altering CDCA and CA. In addition, if this reviewer is understanding the pathway correctly, cbh is catalyzing the deconjugation of conjugated bile acids into free bile acids. If this is the case, would it not be predicted that the adapted strain treated mice show increased free bile acids CDCA and CA rather than less? This deserves additional discussion.\nDiscussion: the statement “the adapted strain led to superior bile acid turnover…and significantly greater protection against disease progression compared to parental strain” is not altogether convincing. It would be helpful to provide some clarity to this comment on bile acid turnover and non-alcoholic steatohepatitis.\nMinor comments:\nSex, age of mice used for the in vivo induction strategy and for the intervention diets should be provided.\nThis reviewer understands that the germ-free mice were employed for the induction strategy. How many mice were used for this part of the study, and from how many of these mice and from what time points were the 60 clones randomly selected from fecal samples collected? This information would be helpful to gauge the scale and scope of the in vivo studies taken."}, "author_response": "Response: We sincerely thank Reviewer #1 for reviewing our manuscript and for their positive feedback on our work. In the following sections, we will address each comment and suggestion made by the reviewer, and the changes in the manuscript are highlighted in yellow.\n1.1 A flowchart is needed to illustrate the process of selecting H22B656.\nResponse to 1.1: Thank you very much for your valuable comments. We agree that a detailed experimental workflow can greatly clarify the study design. In response, we have revised Figure 1a by incorporating a more explicit schematic diagram to illustrate the stepwise screening process of the target Bifidobacterium strain. As shown in the \"In vitro Screening\" panel, we evaluated 12 infant-derived strains based on four key criteria: dissociation ability, enzyme activity, tolerance, and growth performance. A comprehensive multidimensional linear scoring system was then applied to all strains, enabling a robust and quantitative comparison. Through this stringent selection process, Bifidobacterium animalis subsp. lactis H22B656 was ultimately identified as the most promising probiotic candidate for subsequent studies.\nThe specific modifications are as follows:\n1) Revised Figure 1a (line 1152-1164):\nFigure 1. Experimental design and functional screening of infant-derived Bifidobacterium strains. (a) The experimental design comprises three successive phases: In vitro Screening: Twelve infant-derived strains were initially evaluated based on four criteria: dissociation ability, enzyme activity, tolerance, and growth kinetics. This process identified Bifidobacterium animalis subsp. lactis H22B656 as the original candidate strain. Host-mediated In vivo Adaptation: Germ-free mice served as an evolution bioreactor. Following colonization with H22B656, mice were subjected to increasing dietary fat content (30% to 60%) over 7 weeks. Sixty colonies were isolated from feces at different time points. Based on the specific activity ratio, the adapted strain B. animalis subsp. lactis W5S9 was selected. NAFLD Efficacy Validation: The therapeutic potential was evaluated in SPF mice under a high-fat diet (HFAD). Mice were divided into four groups (Control, HFAD + Saline, HFAD + H22B656, and HFAD + W5S9) for an 8-week intervention. Efficacy was assessed via biochemistry, histology, targeted bile acid (BA) metabolomics, and metagenomics to compare the metabolic regulation capabilities of the original and adapted strains.\n1.2 The feces of model animals should be quantitatively tested (RT-qPCR? or other methods...) to reflect the stability of bacterial colonization.\nResponse to 1.2: Thank you for your valuable suggestion. We fully agree that quantitative assessment of bacterial colonization can strengthen the evaluation of strain stability during in vivo intervention. In accordance with this recommendation, we have now incorporated quantitative measurements of Bifidobacterium animalis subsp. lactis H22B656 and its adapted derivative W5S9 in fecal samples. In this study, two experimental components are directly related to bacterial colonization. Below, we describe how colonization was established and quantitatively assessed in each component:\nFirst, in the strain adaptation experiment, to exclude the interference of endogenous microbiota, we used germ-free mice as evolutionary biological reactors to establish a colonization model of Bifidobacterium animalis subsp. lactis H22B656. Each mouse was gavaged with 200 μL of bacterial suspension containing 10⁹ CFU daily for 3 consecutive days, thereby establishing a single-strain colonization model in the germ-free mice.\nSecond, to evaluate the beneficial potential of the adaptive strain Bifidobacterium animalis subsp. lactis W5S9 in modulating host metabolism, we established a diet-induced non-alcoholic fatty liver disease (NAFLD) model in SPF mice. The animals were randomly assigned to four groups (n = 6 per group): control, HFAD, HFAD-H22B656, and HFAD-W5S9. The control and HFAD groups were gavaged daily with 200 μL of sterile saline; the HFAD-H22B656 group received 200 μL of sterile saline containing 10⁹ CFU of the original strain B. animalis subsp. lactis H22B656; and the HFAD-W5S9 group received 200 μL of sterile saline containing 10⁹ CFU of the adaptive strain W5S9. The intervention lasted for 8 weeks. At the end of the experiment, fecal samples were collected for metagenomic sequencing, and taxonomic profiling was performed using MetaPhlAn4 to determine the relative abundance of B. animalis subsp. lactis. Notably, this species was absent from the native intestinal microbiota of the SPF mice prior to intervention, allowing us to specifically track the administered strains. The results showed that both H22B656 and W5S9 established stable colonization in the intestinal tract during the intervention, with comparable colonization efficiency (Fig. 5b in the manuscript).\nFurthermore, based on the reviewer’s suggestion, we additionally quantified the fecal abundance of H22B656 and W5S9 using strain-specific qPCR. The qPCR results were highly consistent with the metagenomic annotations and further confirmed that both the original and adaptive strains achieved stable colonization in the mouse intestine, with similar colonization efficiency (Fig. 5c). We have incorporated this new content into the revised manuscript. The specific revisions are as follows:\n1) Revised Results section (Line 326):\nMetagenomic analysis confirmed that there was no cross‐contamination in the control and HFAD groups, and that oral gavage efficiency was equivalent in the HFAD-H22B656 and HFAD-W5S9 groups, consistent with the fecal qPCR quantification results (Fig. 5b, c; Fig. S2a).\n2) Revised Materials and Methods Section (Line 839-901):\nqPCR Quantification of B. animals in Fecal Samples\nTo quantitatively assess the abundance of the original strain H22B656 and the adapted strain W5S9 in fecal samples, we designed H22B656-specific primers using the NCBI Primer-BLAST tool (H22BF: ACGGCACGATTCATGTGTCT; H22BR: CCAGATTCCATCCTCCTCGC). Fecal samples (0.1 g) were then subjected to DNA extraction using the TIANGEN Bacterial Genomic DNA Extraction Kit, and quantitative real-time PCR (qPCR) was performed with the above primers. qPCR data were analyzed and visualized using GraphPad Prism 10, and statistical significance was evaluated using the Kruskal–Wallis test followed by Dunn’s multiple comparisons test (Fig. 5c).\n3) Revised Figure 5b, c (Line 1256-1260):\n(b, c) Assessment of Bifidobacterium animalis colonization in feces at week 8 via (b) relative abundance analysis (WMS) and (c) qPCR quantification (Log CFU/g). Data are shown as bar plots with individual points (n = 6 mice per group). Kruskal–Wallis test followed by Dunn’s test with multiple comparison correction was used to assess statistical significance;\n1.3 The X-axis of Figure 3c? The Log2(FC)?\nResponse to 1.3: We thank the reviewer for this helpful comment and apologize for the oversight; we have now added the x-axis label (“Log2(FC)”) to Figure 3c.\n1) Revised Figure 3c:\n1.4 Is there any data supporting the impact of H22B656 on the intestinal wall of animals (The level of IL-6? 10? TNF-a?)?\nResponse to 1.4: Thank you for your insightful comment regarding the impact of Bifidobacterium animalis subsp. lactis H22B656 on intestinal barrier integrity and local immune status. In this study, we evaluated the anti-inflammatory efficacy of the original strain H22B656 and the adapted strain W5S9 by quantifying serum levels of key pro-inflammatory cytokines (IL-6 and TNF-α).\nIt is well established that in high-fat diet (HFD)-induced obesity models, systemic inflammation is primarily driven by impaired intestinal barrier function and the subsequent translocation of endotoxins (e.g., lipopolysaccharide, LPS) into the circulation 1. This establishes a direct correlation between elevated serum inflammatory markers and intestinal barrier dysfunction 2, 3. Therefore, although intestinal tissue was not directly sampled in our experiment, the aforementioned pathophysiology supports the use of the serum cytokine profile as a robust proxy for both intestinal barrier integrity and local inflammatory status. Furthermore, elevated levels of IL-6 and TNF-α are well documented to play a pivotal role in disrupting the intestinal mucosal barrier and increasing intestinal permeability 4. Specifically, TNF-α has been shown to promote the degradation of tight junction proteins, such as occludin and members of the claudin family, thereby compromising barrier integrity 2, 3. Concurrently, IL-6 increases epithelial permeability by upregulating the pore-forming protein claudin-2, further impairing the intestinal barrier 4. Therefore, reducing the levels or activity of these pro-inflammatory cytokines is critical for preserving intestinal structural integrity.\nIn our study, supplementation with both the original strain H22B656 and the adapted strain W5S9 significantly reduced serum IL-6 and TNF-α levels compared with the HFAD group (p < 0.05), restoring cytokine levels to those observed in the control group (Fig. 5f). Based on the aforementioned mechanisms, this marked decrease in systemic inflammatory markers strongly suggests that H22B656 mitigates inflammatory stress and thereby preserves intestinal barrier integrity. Although no significant difference was detected between W5S9 and H22B656 in terms of cytokine regulation, both strains exerted significantly stronger anti-inflammatory effects than the HFAD group. The comparable efficacy of W5S9 and H22B656 indicates that the domestication process did not compromise the strain’s functional activity. Collectively, these findings underscore the therapeutic potential of H22B656 in modulating gut-associated immune responses and promoting intestinal health.\nFig.5f Boxplots of serum biomarkers at week 8 (median ± IQR, n = 6). Both H22B656 and W5S9 significantly reduced IL-6, and TNF-α versus HFAD (*p < 0.05; ** p < 0.01; ns., not significant).\nReference\n1. Cani PD, et al. Metabolic endotoxemia initiates obesity and insulin resistance. Diabetes 56, 1761-1772 (2007).\n2. Suzuki T. Regulation of intestinal epithelial permeability by tight junctions. Cell Mol Life Sci 70, 631-659 (2013).\n3. Turner JR. Intestinal mucosal barrier function in health and disease. Nat Rev Immunol 9, 799-809 (2009).\n4. Al-Sadi R, Boivin M, Ma T. Mechanism of cytokine modulation of epithelial tight junction barrier. Front Biosci (Landmark Ed) 14, 2765-2778 (2009).\n1.5 At the final time point, were intestinal segments tissue responsible for metabolizing secondary bile acids collected? And to assess the activation status of relevant targets and the network between bacteria, fermented substances, and the host?\nResponse to 1.5: Thank you for your valuable comments. We would like to begin by clarifying a key concept regarding secondary bile acid metabolism. The biotransformation of primary into secondary bile acids (e.g., deconjugation and 7α-dehydroxylation) is predominantly driven by enzymes secreted by gut microorganisms (such as bile salt hydrolase) within the intestinal lumen, rather than by intestinal epithelial tissues themselves 5, 6. Therefore, to accurately investigate this critical metabolic process, we focused our analyses on colonic contents (fecal samples), which constitute the actual microenvironment where microbial fermentation and secondary bile acid production occur.\nUnfortunately, intestinal tissues were not collected at the final time point of this study for the assessment of local bile acid metabolism or the activation of receptors such as FXR and TGR5. To systematically evaluate the “bacteria–metabolites–host” network, we employed a multi-omics strategy. Shotgun metagenomic sequencing was used to characterize the microbial taxa involved in bile acid metabolism and to annotate key functional genes (e.g., BSH), thereby elucidating the potential mechanisms of bile acid transformation at both the taxonomic and functional levels. Targeted bile acid metabolomics was applied to accurately quantify primary and secondary bile acids. Given that secondary bile acids can influence whole-body metabolic status through the gut–liver axis, we operationally defined “host target activation” as alterations in systemic physiological and biochemical parameters, as well as hepatic histopathology, to reflect the distal regulatory effects mediated via the gut–liver axis.\nIn response, we expanded and refined the correlation analysis and constructed multi-omics association networks integrating microbial taxa, functional genes, bile acid metabolites, and host metabolic indicators in both the original and adapted strain groups (Spearman rank correlation, r > 0.6, p < 0.05). These results are summarized in a new supplementary figure (Fig. S8), which illustrates the expanded multi-omics correlation patterns and provides a more systematic view of their interrelationships. In the HFAD-H22B656 group, Bifidobacterium animalis subsp. lactis H22B656 showed significant correlations with 59 intestinal microbial taxa, indicating robust colonization capacity. However, its regulatory network along the “microbe–metabolite–host” functional axis appeared relatively limited. It was primarily associated with secondary bile acids such as TDCA (taurodeoxycholic acid ， r=-0.6989,p=0.0167) and 7-KLCA (7-ketolithocholic acid, r=-0.6097, p=0.0464), and exhibited significant negative correlations with pro-inflammatory cytokines (TNF-α, r=-0.7485, p=0.0081 and IL-6 r=-0.8080, p=0.0026), suggesting a degree of anti-inflammatory potential. Nonetheless, its correlations with key indicators of hepatic injury (such as AST) were relatively weak, indicating that the original strain had limited capacity to reshape the bile acid profile or ameliorate NAFLD-related metabolic damage. In contrast, the domesticated strain W5S9 exhibited pronounced functional rewiring within its network structure. Most notably, W5S9 showed a strong positive correlation with the BSH functional gene, indicating enhanced bile salt-hydrolyzing potential (BSH, r=0.6681, p=0.0176). Moreover, the bile acid profile of the domesticated-strain group underwent qualitative shifts, with W5S9 exhibiting significant negative correlations with key bile acids such as TCDCA (Taurochenodeoxycholic acid, r=-0.6905, p=0.0129) and TCA (Taurocholic acid, r=-0.6569, p=0.0203), suggesting a stronger modulatory influence on bile acid transformation pathways. At the host level, W5S9 not only retained its inhibitory correlations with inflammatory cytokines (TNF-α , r=-0.6532, p=0.0213 and IL-6, r=-0.7241, p=0.0078), but also showed significant negative correlations with key NAFLD-related pathological indicators (AST, r=-0.6799, p=0.01519 and LDL-C, r=-0.69429, p=0.0123), indicating a markedly enhanced capacity to modulate the bile acid profile and improve host metabolic outcomes.\nIn summary, by performing a systematic correlation analysis that integrates metagenomics, targeted bile acid metabolomics, and systemic host indices, we have constructed a comprehensive multi-omics network. This network holistically elucidates the interaction patterns among the gut microbiota, bile acid metabolism, and host metabolic responses. While we acknowledge the absence of data regarding local receptor activation in intestinal tissues, the current multi-dimensional dataset provides robust evidence to support our key inferences regarding the “Bacteria-Substance-Host” axis. We have explicitly addressed this limitation in the Discussion section and plan to further validate local signaling pathways via transcriptomics and proteomics in future studies.\nReference\n5. Ridlon JM, Harris SC, Bhowmik S, Kang DJ, Hylemon PB. Consequences of bile salt biotransformations by intestinal bacteria. Gut Microbes 7, 22-39 (2016).\n6. Yao L, et al. A selective gut bacterial bile salt hydrolase alters host metabolism. Elife 7, (2018).\nThe specific modifications are as follows:\n1) Revised Discussion section (Line 483-486):\nThrough systematic multi‐omics analysis, we dissected the molecular basis of enhanced bile acid metabolism and ultimately validated the clinical potential of our adapted strains in an animal model to combat NAFLD (Fig. S8).\n2) Revised Figure S8:\nFigure S8. Co-occurrence network analysis revealing the interaction patterns of Bifidobacterium animalis subsp. lactis with the gut ecosystem and host physiological parameters. The networks illustrate the specific correlations between B. animalis subsp. lactis and gut microbiota species, serum biochemical indices, and metabolites in the (a) HFAD-H22B656 and (b) HFAD-W5S9 groups. Each network illustrates Spearman correlation–based associations (|r| > 0.4, p < 0.05) among microbial species (blue nodes), serum biochemical parameters (yellow nodes), bile acid metabolites (green nodes), and microbial functional genes (purple nodes). Positive and negative correlations are indicated by red and blue edges, respectively, with edge transparency corresponding to the magnitude of the correlation.\n1.6 Line 750: The Metaphlan4?\nResponse to 1.6: We thank you for pointing out this issue. We have re-checked the software used in our analyses and confirm that the correct version is MetaPhlAn version 4.0.67. Accordingly, we have corrected the MetaPhlAn version information in the manuscript and updated the reference to the appropriate citation:\n“Blanco-Míguez, A., Beghini, F., Cumbo, F. et al. Extending and improving metagenomic taxonomic profiling with uncharacterized species using MetaPhlAn 4. Nat Biotechnol 41, 1633–1644 (2023).”\nBoth the in-text description and the reference list have been revised to ensure accurate software information and citation.\nReference\n7. Blanco-Miguez A, et al. Extending and improving metagenomic taxonomic profiling with uncharacterized species using MetaPhlAn 4. Nat Biotechnol 41, 1633-1644 (2023).\nThe specific modifications are as follows:\n1) Revised Materials and Methods section (line 885-887):\nTaxonomic profiling of the WMS data was performed using MetaPhlAn (v 4.0.6)7 to identify microbial species and estimate their relative abundances.\n2) Revised Reference section (line 1145-1146):\n67. Blanco-Miguez, A., et al. Extending and improving metagenomic taxonomic profiling with uncharacterized species using MetaPhlAn 4. Nat Biotechnol 41, 1633-1644 (2023).\n1.7 Figure6: The authors only administered H22B656 to germ-free mice. Why is the Shannon diversity index so high? Please provide specific changes in the intestinal bacterial composition.\nResponse to 1.7: Thank you for this important question. We fully understand the reviewer’s concern regarding the unexpectedly high Shannon diversity index shown in Figure 6. We would like to clarify that the data presented in Figure 6 were derived from the efficacy-validation experiment conducted in SPF (specific pathogen-free) mice, rather than from the germ-free (GF) mice used for in vivo strain adaptation. To avoid confusion, we hereby clarify that the current study design comprises two independent animal experiments:\n1. Host-mediated adaptive induction strategy for probiotic strains (using Germ-Free mice)\nGerm-free (GF) mice were employed to establish an “in vivo evolutionary bioreactor” to eliminate interference from any endogenous microbiota. Monocolonization was established via oral gavage of 10⁹ CFU of B. animalis subsp. lactis H22B656 daily for three consecutive days, followed by a 7-week in vivo adaptation process under a stepwise high-fat, high-cholesterol diet. As the mice remained monocolonized throughout this experiment, this phase naturally did not yield data showing elevated Shannon diversity.\n2. Efficacy Evaluation of the Adapted Strain (Using SPF Mice; Source of Data for Figure 6)\nTo evaluate the preventive effect of the adapted probiotic B. animalis subsp. lactis W5S9 on non-alcoholic fatty liver disease (NAFLD), we conducted an experiment using specific pathogen-free (SPF) mice. In this experiment, mice were randomly assigned to four groups (n=6 per group): Control Group: Fed a normal diet and administered 200 μL of sterile saline via oral gavage; HFAD Group: Fed a high-fat diet (HFAD) and administered 200 μL of sterile saline via oral gavage; HFAD-H22B656 Group: Fed a high-fat diet and administered 200 μL of sterile saline containing 109CFU of the original strain B. animalis subsp. lactis H22B656; HFAD-W5S9 Group: Fed a high-fat diet and administered 200 μL of sterile saline containing 109 CFU of the adapted strain B. animalis subsp. lactis W5S9. Following an 8-week intervention, fecal samples were collected for metagenomic sequencing, while serum biochemical indices and liver pathological changes were assessed. Notably, a significant increase in the Shannon diversity index was primarily observed in the HFAD-H22B656 and HFAD-W5S9 groups, indicating that probiotic colonization significantly improved gut microbiota diversity against the background of a high-fat diet (Figure 6a in the manuscript).\nTo further elucidate specific alterations in the bacterial communities, we conducted a detailed taxonomic analysis of the gut microbiota at both the genus and species levels. These results are presented in the newly added Supplementary Figure S3 In terms of overall composition (Fig. S3b), the HFAD group exhibited a distinct shift in microbial structure compared to the Control group. However, intervention with strains H22B656 and W5S9 effectively remodeled the gut microbiota, restoring the compositional profile to a state resembling that of the Control group, particularly regarding the relative abundances of Lachnospiraceae and Muribaculaceae. Further analysis at the species level revealed that the HFAD phenotype was characterized by a depletion of potentially beneficial taxa, such as members of Muribaculaceae (e.g., Muribaculaceae bacterium Isolate 104, Muribaculaceae intestinale) and Paramuribaculum intestinale, concomitant with an enrichment of certain unclassified taxa (e.g., GGB28695, GGB27768). Notably, treatment with both H22B656 and W5S9 successfully reversed these dysbiotic patterns, restoring the abundance of species associated with Muribaculaceae and Bacteroidaceae to levels comparable to the Control group (Fig. S3c, d). We subsequently compared the differential taxa between the two probiotic intervention groups, and differential abundance analysis revealed no significant differences in the abundance of key species, including Bifidobacterium animalis (q-value > 0.05; Fig. S3e). This suggests that both strains exert similar beneficial regulatory effects on gut microbiota composition. Finally, to prevent potential misinterpretation and enhance clarity, we have included a comprehensive experimental schematic in Figure 1, which details the workflow of the study design.\n1) Revised Results section (line 382-383):\nMoreover, among 423 species shared by the two groups, no significant differences in relative abundance were detected (multiple Mann–Whitney tests, q > 0.05; Fig. S3e), indicating that the adapted strain exerts a broadly similar impact on community structure as the original strain, and that W5S9’s functional enhancement did not come at the expense of major shifts in gut microbiota (Fig. S3b, c, d).\n2) Revised Figures section:\nFigure 1(line 1152-1164):\n(a) The experimental design comprises three successive phases: In vitro Screening: Twelve infant-derived strains were initially evaluated based on four criteria: dissociation ability, enzyme activity, tolerance, and growth kinetics. This process identified Bifidobacterium animalis subsp. lactis H22B656 as the original candidate strain. Host-mediated In vivo Adaptation: Germ-free mice served as an evolution bioreactor. Following colonization with H22B656, mice were subjected to increasing dietary fat content (30% to 60%) over 7 weeks. Sixty colonies were isolated from feces at different time points. Based on the specific activity ratio, the adapted strain B. animalis subsp. lactis W5S9 was selected. NAFLD Efficacy Validation: The therapeutic potential was evaluated in SPF mice under a high-fat diet (HFAD). Mice were divided into four groups (Control, HFAD + Saline, HFAD + H22B656, and HFAD + W5S9) for an 8-week intervention. Efficacy was assessed via biochemistry, histology, targeted bile acid (BA) metabolomics, and metagenomics to compare the metabolic regulation capabilities of the original and adapted strains.\nFigure S3\nFigure S3. Effects of B. animalis subsp. lactis H22B656 and B. animalis subsp. lactis W5S9 interventions on intestinal microbiota composition. (a) The Venn diagram depicts the overlap of SNVs called by whole‐metagenome shotgun (WMS) sequencing in the HFAD-W5S9 group and by whole‐genome sequencing (WGS) of the adapted W5S9 strain. Two key SNVs were annotated in the HFAD-W5S9 samples but were not detected in the HFAD-H22B656 group, indicating that no cross-contamination occurred between the gavaged strains. (b) Stacked area chart showing the relative abundance of bacterial taxa at the family and genus levels across the HFAD, HFAD-H22B656, HFAD-W5S9, and Control groups. (c) Venn diagram illustrating the number of shared and unique species between the HFAD-H22B656 and HFAD-W5S9 groups. (d) Heatmap showing the relative abundance of specific bacterial species across all groups. The color scale represents the normalized abundance (Z-score), with red indicating higher abundance and blue indicating lower abundance. (e) Volcano plot comparing the differential abundance of gut microbiota between the HFAD-H22B656 and HFAD-W5S9 groups. The red dot highlights Bifidobacterium animalis, indicating no significant difference (q-value > 0.05) between the two probiotic treatments.\n1.8 The heatmap of 7d-g were scale of row or col?\nResponse to 1.8: We thank the reviewer for the question. In Figures 6 d–g, the heatmaps were scaled by row, i.e., each metabolite was normalized across different groups. This allows for a clearer comparison of the relative changes of each metabolite between groups, independent of the absolute differences across metabolites. The log2FC values presented alongside further indicate the direction and magnitude of changes.\n1.9 The specific software parameters should be provided.\nResponse to 1.9: Thank you for this important suggestion. In the revised manuscript, we have systematically added software versions and key analysis parameters in the Materials and Methods and relevant figure legends. Specifically:\n1) Revised Materials and Methods Section：\nLine 712-716：\nDifferential abundance analysis between the original and adapted strains was performed in R (v4.4.1) using the “DESeq2” (v1.44.0) package (n = 3 per strain, Benjamini–Hochberg FDR correction; FDR < 0.001 and |log₂FC| > 0.4), and the resulting differential features were visualized as volcano plots using OmicShare (https://www.omicshare.com/; Fig. 3c, d).\nLine 804-805：\nStatistical analyses were performed using R v4.4.1 software with the “rstatix” (v0.7.2) and “ggpubr” (v0.6.0) packages.\nLine 885-887：\nTaxonomic profiling of the WMS data was performed using MetaPhlAn (v 4.0.6)7 to identify microbial species and estimate their relative abundances.\nLine 887-891：\nTo specifically distinguish the direct enzymatic contribution of the supplemented strains from potential regulatory effects on the resident microbiota, we performed functional profiling of the metagenomic data using HUMAnN 3 (v3.9) with the uniref90_ko regrouping option to quantify the abundance of bile salt hydrolase (EC 3.5.1.24)–encoding genes (Fig. S6).\nLine 923-932：\nFor targeted fecal bile acid metabolomics, orthogonal partial least squares discriminant analysis (OPLS-DA) was performed using the R (v4.4.1) package “ropls” (v1.36.0) on four groups of mice (n = 6 per group), and differential metabolites were defined as those with a variable importance in projection (VIP) > 1 and log₂ fold change (log₂FC) > 1 between the indicated groups. Score plots and heatmaps were generated with “ggplot2” (v3.5.2) and “pheatmap” (v1.0.12), respectively (Fig. 6d-g). Intersection analysis of differential metabolites was conducted with the Evenn online tool (https://www.bic.ac.cn/test/venn/#/) and visualized as a Venn diagram (Fig. 6h). KEGG pathway enrichment results were displayed as bubble plots using “ggplot2” (v3.5.2) (Fig. 6i), and quantitative data for key differential metabolites were plotted as bar graphs in GraphPad Prism v10 (Fig. 6k).\nLine 936-942：\nTo integrate twelve candidate probiotic strains across multiple performance metrics, we applied a percentile‐ranking linear scoring method and visualized the results as radar plots using the “ggradar” (v0.2) package in R (Fig. 1a). Calcium‐solubilization halo diameters and BSH specific activities for the candidate strains were plotted as bar graphs in GraphPad Prism v10 (Fig. 1b, c). Survival rates under 0.3% bile acid and pH 2 conditions were visualized as heatmaps using the “pheatmap” (v1.0.12) package in R (Fig. 1d), and growth curves were generated with “ggplot2” (v3.5.2) (Fig. 1e).\nLine 949-952：\nStatistical analyses and visualizations, including line and box plots with significance annotations, were generated in R using “ggplot2” (v3.5.2) and “ggsignif” (v0.6.4), with between-group differences assessed by two-tailed Mann–Whitney tests (Fig. 5d, e, f).\nLine 953-959：\nTo accurately describe the regulatory effects of the adapted and original strains on intestinal microbiome composition and bile acid metabolism in NAFLD mice, we used the R package “vegan” (v2.6-10) to calculate α-diversity (Shannon index and species richness) and perform principal coordinates analysis (PCoA) on species-level taxonomic profiles from WMS data (four groups, n = 6 mice per group). Diversity indices and ordination results were visualized as boxplots and scatter plots using “ggplot2” (v3.5.2), and between-group differences in α-diversity indices were assessed by two-tailed Mann–Whitney tests (Fig. 6a–c).\n\nResponse: We sincerely thank Reviewer #1 for reviewing our manuscript and for their positive feedback on our work. In the following sections, we will address each comment and suggestion made by the reviewer, and the changes in the manuscript are highlighted in green.\nMajor Concerns:\n2.1 Lack of mechanistic causality between mutations and phenotype. While the two SNVs (rs1846941 and rs402499) are associated with enhanced bile acid metabolism, direct evidence of causality is limited. Perform complementation or knockout experiments in the parental strain to confirm that these mutations alone are sufficient to confer the enhanced phenotype.\nResponse to 2.1: Thank you for your thorough review and valuable feedback on our study. Regarding your concern about the “lack of mechanistic causality between the mutations and phenotype,” we understand your concern and have taken additional experimental steps to validate the mechanism of the mutations. To investigate the causal impact of these mutations on the phenotype, we performed gene knockout and complementation experiments in both the parental strain Bifidobacterium animalis subsp. lactis H22B656 and the adapted strain Bifidobacterium animalis subsp. lactis W5S9. These experiments were designed to directly verify whether these two mutations are sufficient to confer the enhanced bile acid metabolism phenotype independently.\nFor rs1846941, we knocked out the cbh promoter region where the mutation is located and then assessed the expression of cbh under three different bile salt conditions (0.3% deoxycholic acid (DCA), cholic acid (CA) and chenodeoxycholic acid (CDCA) to explore whether this promoter region regulates gene expression. The results showed that under all three bile salt conditions (DCA, CA, and CDCA), the expression of cbh in the wild-type adapted strain was significantly higher than in the original wild-type strain, suggesting that this non-coding mutation enhances cbh transcription under bile acid stress (fold changes of 1.65, 1.51, and 1.49 for DCA, CA, and CDCA, respectively; One-way ANOVA with Tukey’s HSD test, p < 0.01; Fig. 4d). More importantly, after knocking out the promoter sequence, cbh expression was abolished, further confirming that the upregulation of cbh in the adapted strain is indeed attributed to this non-coding mutation.\nAdditionally, we performed a resazurin assay to evaluate the hydrolytic activity of bile salts in different strains to assess the impact of the promoter sequence on bile salt hydrolysis (Fig. 4e). The results showed that for glycine-conjugated bile salts, the wild-type adapted strain exhibited significantly higher hydrolytic activity than the wild-type original strain. After knocking out the mutation site, both strains lost the ability to hydrolyze glycine-conjugated bile salts. For taurine-conjugated bile salts, the hydrolytic activity of the adapted strain remained higher than that of the original strain, regardless of whether the strain was wild-type or knockout. These findings demonstrate that the non-coding mutation in the cbh promoter region enhances bile salt hydrolysis, particularly in the hydrolysis of glycine- and taurine-conjugated bile salts, by regulating the expression of the cbh gene.\nFor rs402499, we similarly performed a knockout of the mdr gene at the mutation site and examined mdr expression under three different bile salt conditions (0.3% deoxycholic acid (DCA), cholic acid (CA) and chenodeoxycholic acid (CDCA)) to verify whether this mutation directly affects gene expression levels. Under DCA, CA, and CDCA treatment conditions, the mdr expression in the wild-type adapted strain was significantly higher than in the parental wild-type strain, suggesting that this locus significantly enhances mdr transcription in response to bile salt stress (fold changes of 1.74, 2.13, and 1.62 for DCA, CA, and CDCA, respectively; One-way ANOVA with Tukey’s HSD post-hoc test, p < 0.01; Fig. 4f). Further analysis showed that knockout of the mdr gene led to a marked reduction in its expression, indicating that rs402499 directly contributes to the phenotype by regulating mdr transcription. Additionally, we performed complementation experiments to verify whether the restoration of the mdr gene could recover the phenotype. Under CA, DCA and CDCA treatment conditions, the mdr expression in the complemented adapted strain was significantly higher than that in the wild-type complemented strain (p < 0.01; Fig. 4f). These findings suggest that the rs402499 (N192D) mutation drives adaptive phenotypes by increasing mdr expression under bile acid stress.\nMdr, as a transporter protein, can efflux free bile acids, thereby alleviating the toxicity of bile salt molecules to the strain and improving survival 8. To further validate the causal relationship between mdr gene expression and phenotype, we measured the viable bacterial counts of different strains under bile salt stress (p < 0.01; Fig. 4g). The experimental results showed that, under all three bile salt conditions, the colony counts of the wild-type adapted strain were significantly higher than those of the parental wild-type strain. Notably, under CA and CDCA treatment, the mdr knockout strains exhibited a marked growth disadvantage, whereas the differences between knockout strains were less pronounced under DCA conditions (ns). In the complementation experiment, the complemented adapted strain showed significantly higher colony counts than the wild-type parental strain, further supporting the functional role of the mdr gene under bile salt stress. Overall, the mutation in the mdr gene (rs402499) significantly impacted bile acid metabolism and the growth ability of the strain, further confirming the causal relationship between the mutation and the phenotype.\nIn conclusion, gene knockout and complementation experiments provided functional validation of the causal relationship between the non-coding mutations rs1846941 and rs402499 and the associated phenotypic changes. We found that these mutations not only caused significant changes in gene expression but also directly influenced bile salt metabolism and the survival ability of the strains under bile salt stress. Specifically, the non-coding mutation in the cbh promoter region enhanced bile salt hydrolysis activity, while the mutation in the mdr gene improved bile salt transport and strain tolerance. The supplementary content has been added to the manuscript, with the specific modifications outlined as follows:\nReference\n10. Gueimonde M, Garrigues C, van Sinderen D, de los Reyes-Gavilan CG, Margolles A. Bile-inducible efflux transporter from Bifidobacterium longum NCC2705, conferring bile resistance. Appl Environ Microbiol 75, 3153-3160 (2009).\n1) Revised Results section (Line 249-298):\nTo directly verify its regulatory effect within the genomic context, we knocked out the cbh promoter region containing rs1846941 and measured cbh expression levels under three different bile salt conditions (0.3% DCA, CA, and CDCA). The results showed that under these three different bile salt conditions, the cbh expression in the wild-type adapted strain was significantly higher than in the parental wild-type strain, suggesting that this non-coding mutation enhances cbh transcription under bile acid stress (fold changes of 1.65, 1.51, and 1.49 for DCA, CA, and CDCA, respectively; One-way ANOVA with Tukey’s HSD test, p < 0.01; Fig. 4d). Importantly, after knockout of the promoter sequence, cbh expression was abolished, further confirming that the upregulation of cbh in the adapted strain is indeed driven by this non-coding mutation.Moreover, assessment of bile salt hydrolysis using the ninhydrin colorimetric assay revealed that the adapted strain exhibited significantly higher hydrolytic activity toward glycine-conjugated bile salts compared to the parental strain (fold change = 1.25; one-way ANOVA with Tukey’s HSD test, p < 0.01; Fig. 4e). For taurine-conjugated bile salts, hydrolytic activity was also higher in the adapted strain (fold change = 1.16). Notably, knockout of the promoter region significantly reduced hydrolytic activity toward conjugated bile salts in both strains. Together, these findings suggest that the promoter region containing rs1846941 regulates cbh gene expression and plays a critical role in enhancing the hydrolysis of conjugated bile acids.\nFor rs402499, we similarly performed a knockout of the mdr gene, which carries the mutation, and assessed mdr expression under three different bile salt conditions (0.3% DCA; CA; and CDCA) to determine whether the mutation directly affects gene expression. Under these three bile acid treatment conditions, the expression of mdr in the adapted strain was significantly higher compared to the parental wild-type strain, suggesting that this mutation enhances mdr transcription in response to bile acid stress (fold changes of 1.74, 2.13, and 1.62 for DCA, CA, and CDCA, respectively; One-way ANOVA with Tukey’s HSD post-hoc test, p < 0.01; Fig. 4f). Further analysis revealed that knockout of the mdr gene resulted in a significant reduction in mdr expression, confirming that the mutation directly impacts the phenotype through the regulation of mdr expression. Consistently, heterologous expression assays confirmed that the N192D mutation enhances mdr expression at both transcriptional and translational levels under bile salt stress (Fig. S2a, b). Additionally, we performed complementation experiments to assess whether the restoration of mdr expression could recover the phenotype. Under CA, DCA and CDCA treatment conditions, the complemented adapted strain exhibited significantly higher mdr expression than the complemented wild-type strain (One-way ANOVA with Tukey’s HSD post-hoc test, p < 0.01; Fig. 4f). These findings suggest that the rs402499 (N192D) mutation drives adaptive phenotypes by increasing mdr expression under bile acid stress.\nMdr, as a transporter protein, facilitates the efflux of free bile acids, thereby mitigating bile acid toxicity and improving bacterial survival 8. To further validate the causal relationship between mdr gene expression and the phenotype, we assessed the viable bacterial counts under bile acid stress conditions. The results showed that the colony counts of the adapted wild-type strain were significantly higher than those of the parental wild-type strain under all three bile salt conditions (One-way ANOVA with Tukey’s HSD post-hoc test, p < 0.01; Fig. 4g). Notably, under CA and CDCA treatment, the mdr knockout strains exhibited a marked growth disadvantage, whereas the differences between knockout strains were less pronounced under DCA conditions (ns). In complementation experiments, the complemented adapted strain showed significantly higher colony counts than the complemented wild-type strain, further supporting the functional role of mdr gene expression in bile acid stress. Together, the N192D substitution significantly impacted both the transcriptional and translational levels of the mdr gene under bile stress, conferring a robust expression advantage and increasing the strain’s survival rate.\n2) Revised Materials and Methods Section:\nLine 763-791:\nConstruction of gene knockout and complemented strains\nTo establish a causal relationship between the rs1846941 and rs402499 mutations and the phenotype of interest, we employed a homologous recombination strategy using the temperature-sensitive suicide plasmid pKO403 to generate deletion mutants targeting the mdr and cbh promoter regions in both the parental strain Bifidobacterium animalis subsp. lactis H22B656 and the adapted strain B. animalis subsp. lactis W5S9.Genomic DNA from the respective strains was used as the template to amplify the upstream and downstream homologous arms flanking the mdr and cbh promoter regions by PCR using primers listed in Supplementary Table 2. Amplicons were resolved by agarose gel electrophoresis and purified. The pKO403 plasmid was linearized by double digestion with SalI and BamHI and seamlessly assembled with the purified homologous arms using 2×Ezmax® Universal CloneMix (37 °C, 30 min). Assembly products were transformed into E. coli DH5α competent cells, and recombinant plasmids were identified by colony PCR. Plasmids confirmed by PCR were isolated and electroporated into electrocompetent H22B656 and W5S9 cells (2.1 kV, 200 Ω, 25 μF). Following electroporation, cells were recovered in 1 mL MRS-cys medium under anaerobic conditions at 37 °C for 4 h and plated on MRS-cys agar containing chloramphenicol (10 μg/mL) to select for single-crossover integrants. Chloramphenicol-resistant colonies were cultured in MRS-cys broth supplemented with chloramphenicol (10 μg/mL) at 37 °C for 48 h and subsequently plated on MRS-cys agar with chloramphenicol for anaerobic incubation at 42 °C for 48 h to promote double-crossover recombination. Successful deletion of the target region was confirmed by PCR, and double-crossover deletion mutants were obtained. To eliminate the pKO403 plasmid, verified mutants were passaged for three consecutive generations in antibiotic-free MRS-cys medium.\nFor complementation of the Δmdr mutant, a 1,648-bp fragment containing the full-length mdr gene and its native promoter was PCR-amplified from the genomic DNA of H22B656 or W5S9 using primers listed in Supplementary Table 2 and cloned into SalI/BamHI-digested pKO403. The resulting plasmid was introduced into the Δmdr strain by electroporation, and transformants were selected on MRS-cys agar containing chloramphenicol (10 μg/mL). Successful construction of complemented strains was confirmed by PCR.\nLine 793-811:\nRT-qPCR gene expression analysis\nTo assess the impact of the rs402499 and rs1846941 variant on gene expression, bacterial cultures were grown in MRS-cys medium supplemented with 0.3% (w/v) conjugated sodium deoxycholate, sodium cholate, or sodium chenodeoxycholate. After 12 hours of incubation, bacterial cells were collected for subsequent analysis.Total RNA was extracted with TRIzol reagent and treated with DNase I to remove genomic DNA. cDNA synthesis was performed using the PrimeScript™ RT kit according to the manufacturer’s instructions. The cycle threshold (C )values were normalized to the reference gene rplx to obtain ΔC values (ΔC =C -C ). Relative gene expression was calculated using the comparative C method (2-ΔΔCT), with the Wild-type Bifidobacterium animalis subsp. lactis H22B656 strain serving as the calibrator (relative expression = 1) for each bile salt condition independently.\nStatistical analyses were performed using R software with the “rstatix” and “ggpubr” packages. Data are presented as mean ± standard deviation (SD), with individual data points from biological replicates overlaid (n=3). Comparisons between strains under each bile salt condition were analyzed using one-way analysis of variance (ANOVA) followed by Tukey’s Honest Significant Difference (HSD) post-hoc test for multiple comparisons. Statistical significance was defined as p < 0.05. Significant differences are indicated in the figures as follows: *p < 0.05, **p < 0.01, ***p < 0.001; ns, not significant (Fig. 4d, e).\nLine 814-824:\nBile acid stress survival assay for mdr phenotypic validation\nTo validate the contribution of mdr expression to bile acid tolerance, survival under bile acid stress was assessed by viable colony counting. The parental wild-type strain H22B656, the adapted wild-type strain W5S9, their isogenic mdr knockout mutants, and the corresponding complemented strains were grown overnight under standard culture conditions as described above. Cultures were harvested, washed, adjusted to the same optical density, and equal aliquots were inoculated into fresh medium supplemented with 0.3% (w/v) cholic acid (CA), deoxycholic acid (DCA), or chenodeoxycholic acid (CDCA). After incubation under bile acid stress, cultures were serially diluted and plated on bile acid–free agar. Colony-forming units (CFU) were counted after incubation, and survival was expressed as CFU per mL. Each condition was tested in three independent biological replicates, and data were analyzed by one-way ANOVA followed by Tukey’s HSD post hoc test (Fig. 4g).\n3) Revised Figure4 (Line 1220-1244):\nFigure 4. Functional validation of two key SNVs driving enhanced bile acid metabolism in the adapted strain W5S9. (a, b) Sanger sequencing confirmed the presence of the two critical SNVs in the adapted strain W5S9. (c) Comparison of the mCherry fluorescence intensity of reporter plasmids driven by wild-type and rs1846941 mutation under DCA, CA, and CDCA conditions. Bars show mean ± SEM (n = 5). (d, f) Relative expression levels were determined by RT-qPCR in the presence of CA, GDCA, and DCA. Data were normalized to the Wild-type Bifidobacterium animalis subsp. lactis H22B656 (set to 1) independently for each condition. Bars represent the mean ± SD of biological replicates (n=3), with circles indicating individual data points. Statistical comparisons were performed using one-way ANOVA, followed by Tukey’s honestly significant difference (HSD) post hoc test to correct for multiple comparisons. Significance levels: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001; ns, not significant. (e) BSH activity measured by a ninhydrin colorimetric assay using conjugated bile acids as substrate. Bars represent mean ± SD (n = 3). Statistical significance was evaluated using one-way ANOVA with Tukey’s HSD post-hoc test for multiple comparisons. (g) Log CFU of various strains grown under DCA, CA, and CDCA treatments. 10\nBars represent mean ± SD (n = 3). Statistical significance was determined using one-way ANOVA with Tukey’s HSD post-hoc test for multiple comparisons. (h) Molecular docking–predicted binding affinities (kcal/mol) of WT and N192D mdr proteins with DCA, CA, and CDCA. Lower energy indicates stronger binding. Data are presented as mean ± SEM (n = 10). (i) Structural and docking models of WT and N192D mdr proteins with DCA, CA, and CDCA. Left panels show the AlphaFold3-predicted 3D structures of the WT and N192D variants of mdr. Right panels depict the highest-scoring molecular docking poses with DCA, CA, and CDCA, with hydrogen bonds highlighted by dashed lines and red arrows. Statistical significance for docking energies in (h) was assessed using two-tailed Mann–Whitney tests with Benjamini–Hochberg FDR correction; adjusted q values are reported as *q < 0.05, **q < 0.01, and ***q < 0.001.\n2.2 Insufficient characterization of the adapted strain’s ecological impact. The metagenomic data show no major shifts in community structure between HFAD-H22B656 and HFAD-WSS9 groups, but functional changes in the microbiome (e.g., BSH activity of other taxa) are not assessed. Include functional metagenomics or metatranscriptomics to evaluate whether WSS9 influences the bile acid metabolism of other gut microbes.\nResponse to 2.2: We appreciate your insightful comments regarding functional ecological impact. To specifically distinguish the direct enzymatic contribution of the supplemented strains from potential regulatory effects on the resident microbiota, we performed functional profiling of the metagenomic data using HUMAnN v3.9 to quantify the abundance of bile salt hydrolase (EC 3.5.1.24) encoding genes, and incorporated these new results into Supplementary Figure S6. Our analysis revealed that the total abundance of BSH-encoding genes was markedly higher in the HFAD-W5S9 group compared with the HFAD-H22B656 group, which is consistent with the improved bile acid metabolic phenotype observed in the host fecal targeted bile acid metabolomics (Fig. S6a).\nTo further evaluate whether W5S9 induces functional remodeling among background taxa, we additionally performed a stratified contribution analysis, which partitions BSH activity by microbial origin. Notably, the BSH contribution from the endogenous gut microbiota remained largely unchanged between the two groups (Fig. S6b, c), indicating that the elevated BSH capacity is predominantly driven by the colonization and enzymatic expression of W5S9, rather than by extensive functional rewiring of the resident community. These findings indicate that the marked enhancement of total BSH capacity is predominantly attributable to the direct colonization and enzymatic activity of the adapted W5S9 strain itself, thereby conferring a more efficient bile acid transformation phenotype to the gut ecosystem without inducing substantial dysbiosis or functional disturbance within the endogenous microbial community.\n1) Revised Results Section (Line 400-404):\nConsistent with these metabolomic changes, HUMAnN 3-based functional profiling revealed higher total bile salt hydrolase (EC 3.5.1.24) gene abundance in HFAD-W5S9 than HFAD-H22B656 group, with stratified analysis indicating that this enhancement is primarily driven by W5S9 colonization and enzymatic activity rather than broad remodeling of the resident microbiota (Fig. S6).\n2) Revised Materials and Methods Section (Line 888-892):\nTo specifically distinguish the direct enzymatic contribution of the supplemented strains from potential regulatory effects on the resident microbiota, we performed functional profiling of the metagenomic data using HUMAnN (v3.9) with the uniref90_ko regrouping option to quantify the abundance of bile salt hydrolase (EC 3.5.1.24)–encoding genes (Fig. S6).\n3) Revised Figure S6:\nFigure S6. Metagenomic quantification and source attribution of bile salt hydrolase genes. (a) Total abundance of bsh-encoding genes (EC 3.5.1.24) in the gut microbiome across groups. Data are expressed as log10-transformed Reads Per Kilobase (RPKs). (b) Abundance of bsh genes attributed exclusively to the resident background microbiota (excluding the administered Bifidobacterium animalis). (c) Stratified analysis of BSH sources. The stacked bar chart illustrates that the elevated BSH capacity in the HFAD-W5S9 group is primarily driven by the direct contribution of Bifidobacterium animalis (blue), while the contribution from other taxa (grey) remains comparable to the control and HFAD groups. Data are presented as Mean ± SEM. Statistical comparisons were performed using Mann–Whitney test.\n2.3 Limited assessment of long-term stability and safety. The adapted strain was evaluated in an 8-week intervention, but the genetic stability of the mutations under prolonged colonization or in complex microbiota is not tested. Include a stability assay (e.g., serial passaging in vitro or in vivo without selection pressure) to assess whether the mutations revert or confer unintended traits.\nResponse to 2.3：Thank you for your valuable comments regarding the long-term stability and safety of the strain. To address these concerns, we performed a rigorous long-term serial passage assay followed by whole-genome sequencing (WGS). Specifically, the adapted strain W5S9 were subjected to daily serial passaging for 4 weeks under two conditions: (i) a non-selective medium (MRS broth supplemented with L-cysteine hydrochloride), used to assess stability in the absence of selection pressure; and (ii) a selective medium (MRS broth supplemented with bile salts and L-cysteine hydrochloride), used to mimic gut environmental pressure. Target strains were re-isolated weekly, genomic DNA was extracted, and WGS was performed to monitor potential genetic changes. The whole-genome sequencing data of the isolates have been deposited in the NCBI database under BioProject accession number PRJNA1380636.\nThroughout the 4-week passaging period, the identified key SNVs (rs1846941 and rs402499) in the W5S9 strain were maintained at 100% frequency, regardless of whether the medium was selective or non-selective. No reversion to the wild-type genotype was observed, demonstrating that these mutations possess high genetic stability even in the absence of ecological selection pressure (Figure S4a, b).\nWe also wish to clarify that in therapeutic applications, probiotics are typically administered via continuous or semi-continuous dosing, ensuring that the strain ingested daily remains consistent and stable. We have successfully demonstrated phenotypic persistence in our 8-week in vivo study involving continuous administration (The HFAD-W5S9 group received 200 μL of sterile saline containing 10⁹ CFU of the adapted strain W5S9 daily for eight weeks.). Taken together with the new in vitro stability data, we provide robust evidence that W5S9 retains its beneficial properties without undergoing genetic drift within the relevant intervention timeframe.\n1) Revised Results Section (Line 342-344)：\nIn vivo and in vitro safety and stability assessments further demonstrated that the adapted strain W5S9 retains parental-like safety and genetic stability without detectable off-target effects (Fig. S4, S5).\n2) Revised Materials and Methods Section：\nLine 792-736:\nDetermination of genetic stability of key SNVs\nTo verify the genetic stability of key SNVs in the adapted strain, we conducted a rigorous long-term continuous passage experiment followed by whole genome sequencing (WGS). Specifically, the adapted strain W5S9 were subjected to daily serial passaging for 4 weeks under two conditions: (i) a non-selective medium (MRS broth supplemented with L-cysteine hydrochloride), used to assess stability in the absence of selection pressure; and (ii) a selective medium (MRS broth supplemented with bile salts and L-cysteine hydrochloride), used to mimic gut environmental pressure. Target strains were re-isolated weekly, genomic DNA was extracted, and WGS was performed to monitor potential genetic changes (Figure S4a, b).\nLine 965:\nAll genomic resequencing and metagenomic data generated in this study have been deposited in the NCBI database under the following BioProject accession numbers: PRJNA1281804 (third-generation long-read sequencing), PRJNA1281807 (WGS), PRJNA1380636 (WGS), PRJNA1281849 (RNA-seq), and PRJNA1283298 (WMS).\n3) Revised Figure S4:\nFigure S4. Genetic Stability of rs1846941 and rs402499. Mutation frequency of the key SNVs (a) rs1846941 and (b) rs402499 in the adapted strain over a 4-week period under continuous passage conditions. The mutation frequency was assessed in MRS-cys medium (control) and in MRS-cys medium supplemented with 0.3% deoxycholic acid (DCA) as the bile salt stress (n=5 per week).\n2.4 Incomplete translation to human relevance. The study relies solely on a mouse model. The applicability of this approach to human microbiomes remains speculative. Discuss limitations and propose future steps, such as using humanized gnotobiotic models or organoids, to enhance translational relevance.\nResponse to 2.4: We thank the reviewer for raising this important point regarding translational relevance. While the murine model used here enables controlled mechanistic evaluation of bile acid metabolism and host signaling, we recognize that differences in gut ecology and physiology limit direct extrapolation to humans. We have now expanded the Discussion to clearly acknowledge this limitation and to highlight strategies that enhance the translational bridge to human application.\nSpecifically, although the overall microbiome composition differs between mice and humans, the core molecular pathways involved in bile acid metabolism (e.g., BSH-mediated deconjugation) and host signaling (FXR/TGR5) are evolutionarily conserved across both species 9, 10. This conservation suggests that the principal mechanisms underlying the effects of the adapted strain are likely to be cross-species applicable, thereby providing a biological rationale for its future clinical translation.\nTo strengthen translational relevance, we now provide a concrete future research roadmap involving two complementary human-relevant systems:1. Humanized Gnotobiotic Murine Models: By transplanting fecal microbiota from specific cohorts (e.g., patients with hyperlipidemia) into germ-free mice, we aim to evaluate the colonization dynamics and metabolic regulatory capacity of the adapted strain within a human-derived microbial ecology, thereby more accurately recapitulating the human physiological environment.2. Human Intestinal Organoid Systems: We plan to utilize ex vivo human intestinal organoids to directly assess the strain's impact on epithelial barrier function and relevant signaling pathways. This approach will allow us to circumvent limitations imposed by species-specific differences and enhance the physiological relevance of our translational assessment. Together, these future studies will substantially enhance the human translational foundation of our adapted-strain strategy while preserving the mechanistic strengths of the current work. We believe that the present findings provide an essential proof-of-concept for host-guided probiotic optimization and justify continued clinical-oriented development.\nReference\n9. Joyce SA, et al. Regulation of host weight gain and lipid metabolism by bacterial bile acid modification in the gut. Proc Natl Acad Sci U S A 111, 7421-7426 (2014).\n10. Fleishman JS, Kumar S. Bile acid metabolism and signaling in health and disease: molecular mechanisms and therapeutic targets. Signal Transduct Target Ther 9, 97 (2024).\n1) Revised Discussion Section (Line 508-510)：\nTo bridge the translational gap inherent in murine models, future studies will validate the strain's colonization stability and metabolic efficacy using humanized gnotobiotic mice and human intestinal organoid systems.\n2.5 Statistical and methodological clarifications needed. Some figure captions and methods lack detail (e.g., n-values, statistical tests used in metabolomics, multiple testing correction in RNA-seq). Provide full statistical details in figure legends and methods. Ensure that all comparisons are appropriately corrected for multiple testing.\nResponse: We appreciate your constructive feedback. Regarding the lack of detailed explanation of statistics and methodology, we have comprehensively supplemented and standardized the statistical analysis section in the revised draft. The specific revisions are as follows:\n1) Revised Materials and Methods Section：\nLine 712-716：\nDifferential abundance analysis between the original and adapted strains was performed in R (v4.4.1) using the “DESeq2” (v1.44.0) package (n = 3 per strain, Benjamini–Hochberg FDR correction; FDR < 0.001 and |log₂FC| > 0.4), and the resulting differential features were visualized as volcano plots using OmicShare (https://www.omicshare.com/; Fig. 3c, d).\nLine 924-928：\nFor targeted fecal bile acid metabolomics, orthogonal partial least squares discriminant analysis (OPLS-DA) was performed using the R (v4.4.1) package “ropls” (v1.36.0) on four groups of mice (n = 6 per group), and differential metabolites were defined as those with a variable importance in projection (VIP) > 1 and log₂ fold change (log₂FC) > 1 between the indicated groups.\n2) Revised Figure legend Section：\nLine 1215-1218：\n(f, g) RNA-seq–derived counts per million (CPM) for two key upregulated genes: (f) cbh (Log₂ FC = 0.4802, FDR < 0.001) and (g) mdr (Log₂ FC = 0.5355, FDR < 0.001). Bars show mean ± SEM (n = 3 biological replicates per strain; differential expression significance as in (c, d), DESeq2 with Benjamini–Hochberg FDR correction, FDR < 0.001).\nLine 1220-1244：\nFigure 4. Functional validation of two key SNVs driving enhanced bile acid metabolism in the adapted strain W5S9. (a, b) Sanger sequencing confirmed the presence of the two critical SNVs in the adapted strain W5S9. (c) Comparison of the mCherry fluorescence intensity of reporter plasmids driven by wild-type and rs1846941 mutation under DCA, CA, and CDCA conditions. Bars show mean ± SEM (n = 5). (d, f) Relative expression levels were determined by RT-qPCR in the presence of DCA, CA, and CDCA. Data were normalized to the Wild-type Bifidobacterium animalis subsp. lactis H22B656 (set to 1) independently for each condition. Bars represent the mean ± SD of biological replicates (n=3), with circles indicating individual data points. Statistical comparisons were performed using one-way ANOVA, followed by Tukey’s honestly significant difference (HSD) post hoc test to correct for multiple comparisons. Significance levels: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001; ns, not significant. (e) BSH activity measured by a ninhydrin colorimetric assay using conjugated bile acids as substrate. Bars represent mean ± SD (n = 3). Statistical significance was evaluated using one-way ANOVA with Tukey’s HSD post-hoc test for multiple comparisons. (g) Log CFU of various strains grown under DCA, CA, and CDCA treatments. 10\nBars represent mean ± SD (n = 3). Statistical significance was determined using one-way ANOVA with Tukey’s HSD post-hoc test for multiple comparisons. (h) Molecular docking–predicted binding affinities (kcal/mol) of WT and N192D mdr proteins with DCA, CA, and CDCA. Lower energy indicates stronger binding. Data are presented as mean ± SEM (n = 10). (i) Structural and docking models of WT and N192D mdr proteins with DCA, CA, and CDCA. Left panels show the AlphaFold3-predicted 3D structures of the WT and N192D variants of mdr. Right panels depict the highest-scoring molecular docking poses with DCA, CA, and CDCA, with hydrogen bonds highlighted by dashed lines and red arrows. Statistical significance for docking energies in (h) was assessed using two-tailed Mann–Whitney tests with Benjamini–Hochberg FDR correction; adjusted q values are reported as *q < 0.05, **q < 0.01, and ***q < 0.001.\nLine 1259-1261：\nData are shown as or bar plots with individual points (n = 6 mice per group). Statistical significance for multi-group comparisons was assessed using a Kruskal–Wallis test followed by Dunn’s multiple comparisons test with multiplicity-adjusted p values; (d) Body-weight trajectories (mean ± SEM, n = 6 per group).\nMinor Concerns:\n2.6 Figure 3e: The schematic of bile acid metabolism is helpful but could be better annotated with gene names and pathways.\nResponse to 2.6: We appreciate your constructive suggestion. We have revised the schematic in Figure 3e to provide a more detailed and mechanistically informative annotation of the bile acid metabolic network. Specifically, we have now explicitly labeled the genes encoding the key microbial enzymes involved in primary-to-secondary bile acid conversion (e.g., cbh, baiA, baiB, baiCD, baiE). In addition, we used distinct colors and directional arrows to differentiate the major metabolic routes (deconjugation, 7α-dehydroxylation, and related transformations), and we have updated the figure legend with a more comprehensive description of these pathways and their functional roles.\n1) Revised Figure 3e (line 1204-1215):\n(e) Conjugated primary bile acids are taken up by bacteria and undergo deconjugation by bile salt hydrolase (cbh), releasing taurine or glycine and generating free cholic acid (CA) and chenodeoxycholic acid (CDCA). Free primary bile acids may passively diffuse across the cell membrane or be extruded by efflux transporters (e.g., Mdr, LmrA). Within bacterial cells, oxidoreductive reactions mediated by 7β-hydroxysteroid dehydrogenase (7β-HSDH/hdhA) convert CA and CDCA into intermediates such as UCA and UDCA. Subsequently, multi-step 7α-dehydroxylation catalyzed by the bai gene cluster (baiA/baiB/baiCD/baiE) transforms these intermediates into the secondary bile acids deoxycholic acid (DCA) and lithocholic acid (LCA). Transporters such as BaiG and BaiF mediate the export of secondary bile acids into the extracellular space. Red arrows denote steps involved in primary bile acid processing, whereas blue arrows highlight secondary bile acid biosynthesis.\n2.7 The discussion section could be more concise and focused on the broader implications of the work beyond NAFLD.\nResponse to 2.7: We appreciate the reviewer’s insightful suggestion. In response, we have streamlined the Discussion by reducing repetitive NAFLD-specific descriptions and have revised the concluding paragraph to more clearly emphasize the broader implications of our host-mediated in vivo induction strategy beyond NAFLD. The revised text now explicitly highlights how this approach can be generalized to other probiotic functions and disease indications. The updated paragraph in the Discussion reads as follows:\n1) Revised Discussion section (line 487-505):\nImportantly, this host-mediated in vivo induction strategy is not restricted to NAFLD but can be broadly generalized to other probiotic functions and disease indications. By rationally tuning host background, diet, and environmental stressors, the gut can be used as a programmable evolutionary bioreactor to enrich for strains that are optimized for distinct therapeutic tasks. For instance, strains could be adapted in inflammatory or epithelial-injury models to bolster gut barrier restoration and mucosal immune regulation, thereby potentially benefiting inflammatory bowel disease, post-antibiotic gut dysfunction, or complications following chemotherapy and radiotherapy11. Likewise, when evolved under dietary fiber-enriched or prebiotic-supplemented conditions, probiotics might acquire high short-chain fatty acid–producing capacity and refined signaling properties, making them attractive adjuncts for metabolic, cardiovascular, or neuropsychiatric disorders in which gut–brain and gut–liver axes are involved12, 13. Furthermore, applying selective pressures linked to xenobiotic degradation or drug metabolism could yield strains with enhanced capabilities for therapeutic compound modulation or detoxification, broadening their utility in pharmaco-microbiomics and microbiome-based drug-response optimization14. Beyond single-trait enhancement, combining this platform with multi-strain consortia or personalized dietary regimens may ultimately enable patient-specific, functionally tailored live biotherapeutic products that extend well beyond liver disease and provide a generalizable route to next-generation microbiome therapeutics.\nReference\n11. Barbaro MR, et al. Lactiplantibacillus plantarum (CECT7484 and CECT7485) and Pedioccoccus acidilactici (CECT7483) enhance actin cytoskeleton and CYP1A1 expression restoring epithelial permeability alterations induced by irritable bowel syndrome mediators. Gut Microbes 17, 2452235 (2025).\n12. Ke X, et al. Synbiotic-driven improvement of metabolic disturbances is associated with changes in the gut microbiome in diet-induced obese mice. Mol Metab 22, 96-109 (2019).\n13. Haghikia A, et al. Dietary Fatty Acids Directly Impact Central Nervous System Autoimmunity via the Small Intestine. Immunity 43, 817-829 (2015).\n14. Shen S, Wang J, Ma C, Chen Y, Ding H, Zhang J. Understanding the \"individual drug reaction\" from the perspective of the interaction between probiotics and lovastatin in vitro and in vivo. Microbiome 11, 209 (2023).\n3.Please clarify whether the adapted strain was tested for any off-target effects (e.g., growth kinetics, biofilm production, antibiotic resistance, virulence factors, or other traits).\nResponse: We fully appreciate your concerns regarding the safety and stability of the adapted strain. First, it is important to highlight that the parental strain, Bifidobacterium animalis subsp. lactis, is formally listed in the “List of Microorganisms Allowed for Use in Food” in China and has been granted Qualified Presumption of Safety (QPS) status by the European Union. Consequently, it possesses a well-established history of safe use in food products and dietary supplements. Secondly, to directly determine whether the adaptation process introduced any off-target effects, we performed a comprehensive evaluation of the adapted strain across four key dimensions: growth performance (growth kinetics), colonization potential (biofilm formation), antibiotic susceptibility (minimum inhibitory concentrations, MICs), and virulence-associated safety.\nWith respect to growth performance, we carried out a 48-h in vitro growth-curve assay. The adapted strain W5S9 exhibited growth kinetics (lag phase, exponential phase, and stationary-phase biomass) that were indistinguishable from those of the parental strain, indicating that metabolic adaptation did not impose a detectable fitness cost on bacterial growth (Fig. S5a). With regard to colonization potential, biofilm formation was quantified after 12 h of incubation using a crystal-violet assay. No significant difference in biofilm biomass was observed between W5S9 and the parental strain H22B656, suggesting that the adaptive evolution process did not compromise the strain’s intrinsic adhesion or colonization properties (Fig. S5b).\nTo assess antibiotic resistance, we first performed genomic annotation of the parental strain using the Comprehensive Antibiotic Resistance Database (CARD)15. This analysis identified two well-characterized tetracycline resistance genes, tet(W) and tet(M), both encoding ribosomal protection proteins that confer resistance to tetracycline-class antibiotics such as tetracycline and doxycycline16 (Fig. S5c). We then determined the minimum inhibitory concentrations (MICs) of tetracycline and doxycycline for both the parental and adapted strains, and observed identical susceptibility profiles (Fig. S5d). These results confirm that the in vivo adaptation process did not alter the antibiotic susceptibility of the adapted strain W5S9. In parallel, we evaluated virulence-associated safety. Whole-genome annotation using the Virulence Factor Database (VFDB) confirmed the absence of functional toxin-encoding genes, including hemolysins and cytolysins, in the genome of the parental strain H22B65617.\nTaken together, our systematic evaluation (encompassing physiology, colonization, antibiotic susceptibility, and genomic safety) confirms that the adaptation process enhanced metabolic function without inducing off-target effects or compromising safety. Detailed results are available in the revised Supplementary Figure S5.\nReference\n15. Alcock BP, et al. CARD 2023: expanded curation, support for machine learning, and resistome prediction at the Comprehensive Antibiotic Resistance Database. Nucleic Acids Res 51, D690-D699 (2023).\n16. Warburton PJ, Amodeo N, Roberts AP. Mosaic tetracycline resistance genes encoding ribosomal protection proteins. J Antimicrob Chemother 71, 3333-3339 (2016).\n17. Liu B, Zheng D, Zhou S, Chen L, Yang J. VFDB 2022: a general classification scheme for bacterial virulence factors. Nucleic Acids Res 50, D912-D917 (2022).\n1) Revised Results Section (Line 342-344)：\nIn vivo and in vitro safety and stability assessments further demonstrated that the adapted strain W5S9 retains parental-like safety and genetic stability without detectable off-target effects (Fig. S4, S5).\n2) Revised Figure S5:\nFigure S5. Safety evaluation of the parental B. animalis subsp. lactis H22B656 and the adapted B. animalis subsp. lactis W5S9. (a) Growth curves over 48 h showing the growth kinetics of the two strains. (b) Crystal-violet biofilm formation assay after 12 h of anaerobic incubation; biofilm biomass is expressed as OD590/OD600. (c) Resistome profile based on annotation with the Resistance Gene Identifier (RGI) against the CARD database; only hits with >95% amino-acid identity to reference sequences are shown. orf_id, open reading frame identifier; Identity (%), percentage amino-acid identity to the reference gene. (d) Dose–response growth-inhibition curves of the two strains in the presence of doxycycline and tetracycline.\n\nResponse: We sincerely thank Reviewer #1 for reviewing our manuscript and for their positive feedback on our work. In the following sections, we will address each comment and suggestion made by the reviewer, and the changes in the manuscript are highlighted in blue.\nMain comments:\n3.1 Figure 4d: The authors state that “mdr expression was significantly higher in the mdr-N192D under all bile acid conditions”. While that might be accurate at the 3h time point, it does not fully describe the data shown about mdr expression going down from baseline for both WT and N192D variant following induction with cholic acid CA and chenodeoxycholic acid CDCA. In addition, the followup sentence “Together, these results demonstrate that the N192D substitution upregulates mdr expression” isn’t entirely accurate…could this statement be revised given that the variant upregulates mdr expression to a greater degree after induction specifically with deoxycholic acid induction and not with CA or CDCA? The second half of the statement “potentially by altering local mRNA secondary structure (thus enhancing ribosome binding or translation initiation), optimizing codon usage bias, modulating mRNA folding dynamics or modifying regulatory elements affecting splicing translation” seems overly speculative.\nResponse to 3.1: We sincerely appreciate your valuable comments and the important points you raised regarding Figure 4d. Regarding your observation of the statement “mdr expression was significantly higher in the mdr-N192D variant under all bile acid conditions,” we truly value your insight. Upon further review of the data, we realize that this statement may not fully reflect the changes in mdr expression, particularly under cholic acid (CA) and chenodeoxycholic acid (CDCA) induction conditions, where both the wild-type and N192D variant showed a decrease in mdr expression compared to baseline levels after induction. We have revised the original statement to: “RT-qPCR analysis revealed distinct substrate-dependent patterns: under DCA, mdr-N192D exhibited upregulation relative to the baseline (0 h), whereas mdr-WT levels remained relatively stable. Conversely, exposure to CA or CDCA led to a repression of mdr expression in both variants; however, the mdr-N192D variant displayed a significantly attenuated reduction compared to the wild-type.”\nAdditionally, I would like to clarify that in order to more clearly elucidate the causal relationship between the mutation and phenotype, we have supplemented the study with gene knockout and complementation experiments. To further optimize the manuscript’s logical structure and ensure a more coherent presentation, we have moved the results from the heterologous expression and reporter system experiments to the supplementary figures (Fig. S8a, b). We have shifted the focus to the direct results from the gene knockout and complementation experiments, which we believe better explain the impact of the N192D mutation on mdr expression and its associated phenotype. This modification should enhance the clarity of our main findings while allowing us to concentrate on more direct experiments that confirm the mutation’s effect on mdr expression and its role in bile acid metabolism.\nFurthermore, we acknowledge that the second part of the statement regarding potential mechanisms (e.g., altering local mRNA secondary structure, optimizing codon usage bias, etc.) was speculative. Based on your feedback, we have removed this speculative portion.\n1) Revised Results Section (Line 276-278):\nConsistently, heterologous expression assays confirmed that the N192D mutation enhances mdr expression at both transcriptional and translational levels under bile salt stress (Fig. S2a, b).\n3.2 Line items 272-273: in the results section, it states “N192D mutation stabilizes the protein-substrate complex, thereby facilitating more efficient substrate recognition and transmembrane transport” is declarative without evidence…would it be more accurate to state “…complex, potentially facilitating more efficient substrate recognition…”?\nResponse to 3.2: We sincerely thank you for this precise suggestion regarding the accuracy of our language. We agree that the original statement was overly declarative, as the conclusion is based on computational simulations rather than direct kinetic evidence of transport. We have adopted your suggestion and tempered the language to more objectively reflect the insights provided by the simulation data. The sentence has been revised to: “...thereby potentially facilitating more efficient substrate recognition and transmembrane transport.”\nThe specific modifications are as follows:\n1) Revised Results Section (Line 307):\nThese results suggest that the N192D mutation stabilizes the protein–substrate complex, thereby potentially facilitating more efficient substrate recognition and transmembrane transport.\n3.3 Adapted strain on weight gain: The three parameters chosen are useful outcome measurements to gauge the effect of the adapted strains. However, what could be clarified is the main comparison of interest. In the weight gain figure (Figure 5c), the authors state the HFAD-W5S9 treated mice showed the smallest weight gain from control, followed by the parental strain treated mice. The effect size of the adapted strain is not immediately accessible to the reader. What is the absolute difference in effect?…did it prevent weight gain imposed by a HFD by 30%, 50%, etc. and how does that compare to the parental strain? The results section could be rewritten to clarify this point.\nResponse to 3.3: We sincerely appreciate your insightful comment. We fully acknowledge that quantifying the effect size allows for a more direct and intuitive comparison of the therapeutic potential between the adapted and parental strains, thereby providing clearer evidence of the intervention effects. To address this point, we reanalyzed the body weight data by calculating both the net body weight gain during the intervention period and the percentage reduction in weight gain relative to the high-fat diet (HFAD) model group. The quantitative analysis revealed that the adapted strain W5S9 significantly attenuated HFAD-induced body weight gain, achieving a 33.6% reduction (mean gain: 7.76 g vs. 11.68 g in the model group; p = 0.0022). Moreover, although the parental strain H22B656 exhibited a numerically reduced gain (8.13 g, 30.4% reduction), the difference did not reach statistical significance compared with the model group (p > 0.05). These results collectively indicate that adaptive evolution confers a measurable efficacy advantage to W5S9 over its parental strain in preventing HFAD-induced weight gain.\nThe specific modifications are as follows:\n1) Revised Results Section (Line 335-340):\nTo quantify this therapeutic effect, we analyzed the net weight gain. W5S9 intervention significantly prevented HFAD-induced weight gain by 33.6% (mean gain: 7.76 g vs. 11.68 g in the HFAD group; Mann-Whitney test, p=0. 0022), exhibiting the smallest weight gain among the HFAD-fed groups. In contrast, the original strain showed a less pronounced reduction with higher variability (mean gain: 8.13 g vs. 11.68 g in the HFAD group; inhibition rate=30.4%, p=0. 1320, Fig. 5e).\n2) Revised Figures 5e (Line 1263-1265):\n(e) Bar plot displaying the net increase in body weight (g) across groups over the 8-week intervention period, with individual measurements shown. Values represent mean ± SEM. Group differences were evaluated using the Mann–Whitney U test (*p < 0.05; ** p < 0.01).\n3.4 Adapted strain on liver histology: It would be helpful to perform semi-quantitative scoring for histology and oil red O staining across samples within each group to get a more accurate assessment of the adapted strain’s impact, if any, on the liver. This is critical as non-alcoholic steatohepatitis rests ultimately on histologic diagnosis. Both the parent and adapted strain appeared to reduce AST levels so it is unclear if there are significant differences in their impact on liver pathology.\nResponse to 3.4: We sincerely appreciate your valuable suggestion, especially considering that the diagnosis of non‑alcoholic steatohepatitis (NASH) relies on liver histopathology. To address this critical issue, we performed quantitative analysis of Oil Red O‑stained liver sections from all experimental groups (n = 6) using ImageJ software, calculating the percentage of lipid droplet area relative to the total liver section area. The results have been incorporated as a new panel in the revised Figure 5i.\nThe quantitative results confirmed and extended our initial observations. Compared with the HFAD group, mice in the HFAD‑W5S9 group showed a significantly lower hepatic lipid burden (p= 0.0225, Kruskal–Wallis with Dunn’s test), reducing lipid accumulation to a level close to that of the control group (Figure 5i), with no statistically significant difference in lipid droplet area compared with healthy controls. This indicates that the adapted strain exerts a more pronounced protective effect against hepatic steatosis.\n1) Revised Results Section (Line 363-366):\nQuantitative image analysis of lipid droplet area (n = 6 per group) further reinforced this observation, revealing a significantly greater reduction in hepatic lipid content in the HFAD-W5S9 group compared to the HFAD group (p= 0.0225, Kruskal–Wallis with Dunn’s test), approaching levels observed in the healthy control group (Fig. 5i).\n2) Revised Materials and Methods Section (Line 870-892):\nSemi-quantitative analysis of lipid droplet content was performed on Oil Red O–stained sections using Fiji (v2.17.0). Data visualization and statistical analyses were conducted in GraphPad Prism v10 using a Kruskal–Wallis test followed by Dunn’s multiple comparisons test (n = 6 mice per group, Fig.5i).\n3) Revised Figure 5i (Line 1270-1274):\n(i) Quantification of lipid accumulation in liver sections stained with Oil Red O. Bar graph shows the mean lipid droplet area per field across groups (n = 6 mice per group). Data are expressed as mean ± SEM. Statistical significance was determined using the Kruskal–Wallis test followed by Dunn’s multiple comparisons test (*p < 0.05; ** p < 0.01, *** p < 0.001, **** p < 0.0001).\n3.5 Figure 6 figure legend title “The adapted B. animalis subsp. Lactis W5S9 re-balances intestinal bile-acid metabolism and prevents NAFLD progression” needs to re-worded, as “prevents NAFLD progression” is not shown in this figure.\nResponse to 3.5: We sincerely thank you for this precise correction. We agree that the original title was an overstatement regarding the specific data presented in Figure 6, which primarily illustrates changes in the gut microbiota and the modulation of bile acid metabolism, rather than the direct phenotypic outcomes of NAFLD prevention (which are demonstrated in Figure 5). Accordingly, we have revised the figure legend title to strictly describe the data shown. The title has been changed to: “The adapted B. animalis subsp. lactis W5S9 re-balances intestinal bile-acid metabolism in NAFLD mice.”\n1) Revised figure legend section (line 1277-1278):\nFigure 6. The adapted B. animalis subsp. lactis W5S9 re-balances intestinal bile-acid metabolism in NAFLD mice.\n3.6 Adapted strain on bile acid metabolites: In figure 6K, the authors note that bile acid concentrations in the adapted strain treated mice approached those of the control group more closely than in the parental strain treated mice, “indicating a more robust reestablishment of bile acid homeostasis”. This reviewer finds it interesting that the probiotic parental strain also impacted bile acid metabolism. The data is presented in such a way that it is difficult to gauge the effect of parental strain vs the adapted strain in altering CDCA and CA. In addition, if this reviewer is understanding the pathway correctly, cbh is catalyzing the deconjugation of conjugated bile acids into free bile acids. If this is the case, would it not be predicted that the adapted strain treated mice show increased free bile acids CDCA and CA rather than less? This deserves additional discussion.\nResponse to 3.6: We sincerely appreciate your valuable comments regarding the bile acid metabolism results. In response to your question about the trends in free bile acids (CDCA and CA) shown in Figure 6K, we have revised the figure to improve clarity. Specifically, we removed unnecessary characters and now only display the log₂ fold change (Log₂FC) of the HFAD-H22B656 and HFAD-W5S9 groups relative to the HFAD group(All the metabolites presented in the updated figure were significantly different from the HFAD group (VIP > 1, |Log₂FC| > 1)).We would like to clarify that the previous figure included a misinterpretation of certain values due to overly compact formatting. This issue was not raised by the reviewer but was discovered by our team during subsequent revision and figure checking. We have now carefully re-verified the data and corrected both the figure and the manuscript accordingly.\nFurthermore, to help visualize the differences between the parental and adapted strains more clearly, we have added Supplementary Figure 7, which presents the Log₂FC values in a more intuitive and informative layout.\nHere, we would like to further elaborate on the differential efficacy between the parental and adapted strains. While the parental strain H22B656 exerts beneficial effects on bile acid metabolism, the adapted strain W5S9 demonstrates a substantially stronger capacity to restore bile acid homeostasis. In particular, W5S9 more effectively reduces the HFAD-induced elevations of conjugated bile acids such as TCDCA and TCA, bringing their levels closer to those of the healthy control group (Fig. 6k; Fig. S7). These findings highlight the functional enhancement conferred by host-mediated adaptation and underscore the superior robustness of W5S9 in correcting HFAD-associated metabolic disturbances.\nIn addition, regarding your second question: since cbh catalyzes the deconjugation of conjugated bile acids into free bile acids, one might expect that treatment with the adapted strain would lead to increased levels of free bile acids such as CDCA and CA. You have indeed raised an important theoretical point. According to classical biochemical pathways, if bile salt hydrolase (BSH) catalyzes the deconjugation of conjugated bile acids (e.g., TCDCA → CDCA), it would be intuitive to assume that the resulting free bile acids (CDCA and CA) should accumulate. However, our data instead show decreased levels of these free bile acids. This seemingly contradictory observation reflects the complexity of enterohepatic circulation and systemic regulatory mechanisms that operate in vivo, which cannot be explained by a simple linear consequence of a single enzymatic reaction18, 19.Although BSH activity indeed promotes the conversion of conjugated bile acids into free forms (e.g., TCDCA → CDCA), the ultimate fate of these deconjugated products is determined by the dynamic control of the enterohepatic system. Under physiological conditions, bile acid homeostasis involves not only enzymatic deconjugation, but also fecal excretion, intestinal reabsorption, hepatic de novo synthesis, and microbial secondary transformations—interconnected processes that cannot be fully reproduced in closed in vitro enzymatic assays20.\nThe observed reduction in fecal CDCA and CA levels is primarily attributable to an enhanced “excretion effect.” A key physiological consequence of bile salt hydrolase (BSH) activity is the inhibition of efficient reabsorption of conjugated bile acids. Under normal conditions, conjugated bile acids are actively reabsorbed by the host via specific transport systems such as the apical sodium-dependent bile acid transporter (ASBT) in the ileum21. However, once deconjugated by BSH, the resulting free bile acids (e.g., CDCA and CA) exhibit markedly lower aqueous solubility and are reabsorbed with substantially reduced efficiency9. Consequently, these free bile acids are more prone to be eliminated via fecal excretion rather than accumulating in the serum or liver22. Under this mechanism, BSH activity does not lead to an “accumulation” of free bile acids; rather, by promoting their fecal elimination, it effectively contracts the overall bile acid pool. It is important to note that our samples were derived from fecal-targeted bile acid metabolomics, and therefore reflect the absolute concentrations of terminal excretory products. Thus, the observed reduction of free bile acids (CDCA and CA) in feces compared to the HFAD group should not be interpreted as indicative of diminished BSH activity or insufficient deconjugation. Instead, it reflects a systemic regulatory outcome of restored metabolic homeostasis and improved bile acid clearance. In this context, as the conjugated bile acid precursors are progressively depleted and the total bile acid pool contracts, the decline in fecal CDCA and CA concentrations represents a physiologically expected result.\nTaken together, the reduced levels of CDCA and CA observed in the adapted-strain group are therefore not indicative of insufficient BSH activity, but rather reflect the successful re-establishment of bile acid homeostasis. By coupling enhanced deconjugation-driven excretion with the restoration of physiological feedback inhibition, the adapted strain effectively contracts the HFAD-induced pathological expansion of the bile acid pool, resulting in a bile acid profile that closely resembles that of healthy controls.\nReference\n9. Joyce SA, et al. Regulation of host weight gain and lipid metabolism by bacterial bile acid modification in the gut. Proc Natl Acad Sci U S A 111, 7421-7426 (2014).\n17. Liu B, Zheng D, Zhou S, Chen L, Yang J. VFDB 2022: a general classification scheme for bacterial virulence factors. Nucleic Acids Res 50, D912-D917 (2022).\n18. Chiang JY. Bile acids: regulation of synthesis. J Lipid Res 50, 1955-1966 (2009).\n19. Sayin SI, et al. Gut microbiota regulates bile acid metabolism by reducing the levels of tauro-beta-muricholic acid, a naturally occurring FXR antagonist. Cell Metab 17, 225-235 (2013).\n20. Ramirez-Perez O, Cruz-Ramon V, Chinchilla-Lopez P, Mendez-Sanchez N. The Role of the Gut Microbiota in Bile Acid Metabolism. Ann Hepatol 16 Suppl 1, S21-S26 (2017).\n21. Dawson PA. Role of the intestinal bile acid transporters in bile acid and drug disposition. Handb Exp Pharmacol, 169-203 (2011).\n22. Begley M, Gahan CG, Hill C. The interaction between bacteria and bile. FEMS Microbiol Rev 29, 625-651 (2005).\n1) Revised Results Section (Line 422-428):\nTo comprehensively evaluate the regulatory magnitude of the strains, we analyzed the Log2 fold changes (Log2FC) relative to the HFAD model (Fig. S7). While the parental strain H22B656 elicited a moderate reduction in pathological bile acid accumulation, the adapted strain W5S9 exhibited a significantly more profound restorative effect. Specifically, W5S9 treatment resulted in a substantial decrease in both conjugated bile acids (TCDCA, log₂FC = –1.14; TCA, log₂FC = –1.19; GCDCA, log₂FC = –∞; GCA, log₂FC = –2.68;) and free bile acids (e.g., CDCA, log₂FC = –2.09; CA, log₂FC = –2. 02) in the intestinal lumen (Fig. 6k).\n2) Revised Figures Section:\nFigure 6k (Line 1296-1299):\n(k) Concentration of the six key bile acids highlighted in (j). Log₂ fold-change vs HFAD is annotated above each panel (All bile acid metabolites showed significant differences compared to the HFAD group (VIP > 1, |Log₂FC| > 1), mean ± SEM, n = 6).\nFigure S7\nFigure S7. Targeted bile acid metabolomics showing log2 fold change in fecal bile acid concentrations relative to the HFAD group. Bar plots represent the mean ± SEM for each bile acid metabolite, with individual dots indicating biological replicates.\n3.7 Discussion: the statement “the adapted strain led to superior bile acid turnover…and significantly greater protection against disease progression compared to parental strain” is not altogether convincing. It would be helpful to provide some clarity to this comment on bile acid turnover and non-alcoholic steatohepatitis.\nResponse to 3.7: We thank the reviewer for this helpful comment. We agree that the original wording was too general and did not sufficiently clarify how bile acid “turnover” was defined or how the model relates to non-alcoholic steatohepatitis (NASH). To address this, we have revised the sentence in the Discussion to (i) specify that the diet-induced NAFLD model recapitulates key features of NAFLD and (ii) replace the vague term “superior bile acid turnover” with a more explicit description of bile-acid homeostasis and the associated biochemical and histological changes. The revised sentence now reads:\n1) Revised Discussion Section (Line 444-450):\nSubsequently, in a diet-induced NAFLD mouse model that recapitulates key features of non-alcoholic steatohepatitis, dietary supplementation with the adapted strain more effectively restored bile-acid homeostasis, as indicated by a stronger reduction and normalization of both conjugated (TCDCA, GCDCA, TCA and GCA) and free primary bile acids (CDCA, CA), and was accompanied by lower inflammatory markers and more pronounced attenuation of hepatic steatosis and liver injury than observed with the parental strain.\nMinor comments:\n3.8 Sex, age of mice used for the in vivo induction strategy and for the intervention diets should be provided.\nResponse to 3.8: We apologize for not explicitly stating the demographic details of the experimental animals in the original manuscript. We have now revised the “Materials and Methods” section to strictly define the sex, age, and strain of the mice used for both the in vivo adaptation phase and the dietary intervention phase.\nFor the in vivo adaptation strategy: We have specified that male GM C57BL/6J mice, aged 7 weeks, were used as the evolution bioreactor.\nFor the intervention diets (NAFLD model): We have clarified that male SPF C57BL/6J mice, aged 7 weeks, were employed for the probiotic intervention study.\n1) Revised Materials and Methods section:\nLine 590:\nMale GM C57BL/6J mice (7 weeks old; n = 4) were obtained from Shenzhen Jingtai Biotech (Guangdong, China) and maintained in sterile isolators on irradiated standard chow for 7 days.\nLine 847-850:\nDiet-induced NAFLD mouse model and probiotic gavage\nTo assess the beneficial potential of the adapted B. animalis subsp. lactis W5S9 in modulating host metabolism, we employed a diet‐induced NAFLD mouse model in male SPF C57BL/6J mice (7 weeks old), based on the established link between bile acid dysregulation and NAFLD progression.\n3.9 This reviewer understands that the germ-free mice were employed for the induction strategy. How many mice were used for this part of the study, and from how many of these mice and from what time points were the 60 clones randomly selected from fecal samples collected? This information would be helpful to gauge the scale and scope of the in vivo studies taken.\nResponse to 3.9：We verify that 4 GM C57BL/6J mice were used as the in vivo evolution bioreactor for this induction strategy. These mice were co-housed in a sterile isolator to maintain a uniform microbial environment. Regarding the selection of the 60 clones: To capture the evolutionary trajectory of the strain under the stepwise increasing selection pressure, fresh fecal samples were collected from all 4 mice at four specific time points: the end of Weeks 1, 3, 5, and 7. At each time point, fecal samples were homogenized, and colonies were isolated on MRS-cys agar. We randomly selected 15 isolates from each time point (totaling 60 isolates: 15 clones × 4 time points) for downstream BSH activity screening. This strategy allowed us to assess the dynamic enhancement of metabolic capacity throughout the induction period.\n1) Revised Materials and Methods section:\nLine 602-606:\n“...To ensure a representative evaluation of the adaptation process, 15 isolates were randomly selected from the plates of each sampling time point (Weeks 1, 3, 5, and 7), resulting in a total library of 60 induced isolates. These isolates were inoculated into MRS-cys liquid medium, labeled sequentially, and stored at –80 °C for subsequent analyses...”\nReviewer #3 (Remarks to the Author):\nI have no further comments. The authors have adequately addressed this reviewer's prior comments.\nThere is a typo in line 138 that refers to Figure 1a. shouth that be Figure S1b?\nResponse to Reviewer #3:\nWe thank the reviewer for their positive comments. We have corrected the typo in line 138, and the figure reference now correctly refers to Figure S1b.\nLine136-139：\n“To comprehensively evaluate probiotic potential in bile acid metabolism, we integrated these five key indicators from four dimensions and conducted a percentile‐based ranking using a multidimensional linear summation score for each strain. The resulting composite scores were visualized as radar plots (Supplementary Fig.1b).”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Congratulations on your valuable study.", "reviewer_2": "The authors have responded comprehensively to all comments. The revised manuscript presents a rigorous, innovative, and well-rounded study that advances the field of probiotic engineering. The experimental additions have substantially strengthened the causal and ecological conclusions, and the expanded discussion thoughtfully contextualizes the work within current translational challenges. I therefore recommend acceptance.", "reviewer_3": "I have no further comments. The authors have adequately addressed this reviewer's prior comments.\nThere is a typo in line 138 that refers to Figure 1a. shouth that be Figure S1b?"}, "author_response": ""}]

Decision letter