Five- to six-week-old male Wistar rats were obtained from Jinan Pengyue Experimental Animal Breeding Co., Ltd. (Jinan, China; certificate no. 370726230100193163). The rats were acclimatized in a specific pathogen-free animal room that was maintained at 22 ± 2°C with 50–70% humidity on a 12-h light/dark cycle and provided water and food.

After acclimatization for 1 week, the eight rats in the control (Con) group were fed standard chow for 8 weeks. In addition, 32 rats were fed a high-fat diet (HFD) for 8 weeks. NAFLD was confirmed via liver pathology in three rats randomly selected from the HFD group. On the basis of their relatively high body weights, 24 rats were randomly divided into three groups (n = 8/group). The Wistar rats in the Con group (n = 8) and the NAFLD model (Mod) group received distilled water via oral gavage. The OCA group was administered 100 mg/kg/d OCA (Shanghai Yuanye Bio-Technology Co., Ltd., China) via gavage (Xin et al., 2021), and the XFD group was administered 18 g/kg/d XFD (First Affiliated Hospital of Anhui University of Chinese Medicine, China) via gavage. This rat XFD dosage is six times greater than the human dosage (Zhao and Fang, 2024). The appropriate treatment was administered once a day, and the treatment course lasted 12 weeks.

XFD is composed of Bambusae caulis in Taeniam (Anhui Puren Traditional Chinese Medicine Pieces Co., Ltd., 10 g), Atractylodis Rhizoma (Bozhou Huqiao Pharmaceutical Co., Ltd., 15 g), Bupleuri Radix (Anhui Xiehecheng Pharmaceutical Co., Ltd., 10 g), Aurantii Fructus (Anhui Puren Traditional Chinese Medicine Pieces Co., Ltd., 15 g), Scutellariae Radix (Anhui Xiehecheng Pharmaceutical Co., Ltd., 10 g), Curcumae Radix (Anhui Puren Traditional Chinese Medicine Pieces Co., Ltd., 15 g), Corydalis Rhizoma (Anhui Boyao Qiancao National Pharmaceutical Co., Ltd., 15 g), Paeoniae Radix Alba (Anhui Xiehecheng Pharmaceutical Co., Ltd., 20 g), Taraxaci Herba (Anhui Xiehecheng Pharmaceutical Co., Ltd., 20 g), Crataegi Fructus (Bozhou Huqiao Pharmaceutical Co., Ltd., 15 g), Setaria italica (L.) Beauv (Bozhou Huqiao Pharmaceutical Co., Ltd., 20 g), and Plantaginis Herba (Bozhou Huqiao Pharmaceutical Co., Ltd., 15 g). XFD was soaked in 1,800 mL of water. The XFD was boiled at high heat and then simmered at low heat for 30 min. The liquid medicine was collected. XFD was subsequently boiled for a second time as previously described. The liquid mixture was mixed and concentrated to 50 mL. The XFD concentration was 3.6 g/mL.

Body weights were measured via a weight scale. After 12 h of fasting, the fasting blood glucose (FBG) level was measured via the tail vein with a Roche ACCU-CHEK Performa glucometer (Basel, Switzerland). The rats were subsequently anesthetized via intraperitoneal injection of 30 mg/kg pentobarbital sodium (Merck, United States). Blood samples were collected from the abdominal aorta and centrifuged at 4,000 revolutions per minute (rpm) for 10 min. Furthermore, the entire liver tissue samples were weighed, and the liver lobe and white fat from the same area were preserved. Feces were removed by dissecting the ileum. The serum, liver, and feces were stored at −80°C.

The liver was weighed, ground, and centrifuged, after which the supernatant was collected. The levels of total cholesterol (TC), triglyceride (TG), alanine aminotransferase (ALT), and aspartate aminotransferase (AST) in the serum and the levels of TG and TC in the liver were determined using a Chemray 800 fully automatic biochemical analyzer (Rayto Life and Analytical Sciences Co., Ltd., China). The fasting insulin (FIns; American Laboratory Products Company, United States) and bile salt hydrolase (BSH; Jiangsu Medical Industrial Co., Ltd., China) concentrations were analyzed via enzyme-linked immunosorbent assay (ELISA) kits according to the manufacturers’ instructions. The homeostasis model assessment–insulin resistance (HOMA–IR) score was calculated as FIns (μIU/mL) × FBG (mmol/L) ÷ 22.5.

The liver and white fat portions were fixed in a tube containing a 4% paraformaldehyde solution for 24 h. After dehydration, the samples were embedded in paraffin wax, sliced into 4 μm sections, and stained with hematoxylin and eosin (HE).

The liver and white fat tissues were embedded in optimal cutting temperature compound to prepare frozen sections. The samples were then cut and immersed in oil red O staining solution for 8–10 min before removal and subsequent immersion in each of two cylinders of 60% isopropanol for 3 s and 5 s for differentiation. The slices were then soaked in each of 2 tanks of pure water for 10 s, sliced, immersed in hematoxylin counterstain for 3–5 min, and then washed in each of 3 tanks of pure water for 5, 10, and 30 s. The samples were subsequently differentiated for 2–8 s, washed in each of two tanks of distilled water for 10 s, returned to the staining solution for 1 s, and then immersed in each of two tanks of tap water for 5 s and 10 s. Finally, the slices were sealed with a glycerol gelatin sealing agent. Aipathwell (Servicebio) software was used to identify the positive and tissue areas within the samples. The ratio of oil red O-stained area was calculated as the stained area/tissue area. Images were captured under an Eclipse Ci-L microscope (Nikon Corporation, Japan).

Total genomic DNA was extracted from the fecal samples using an E.Z.N.A.^® Soil DNA Kit (Omega Bio-tek, United States). The extracted DNA was randomly fragmented to a size of approximately 400 bp using a Covaris M220 instrument (Gene Company Limited, China). A paired-end library was constructed with NEXTFLEX Rapid DNA-Seq (Bios Scientific, United States). Adapters containing the complements of the sequencing primer hybridization sites were ligated to the blunt ends of the fragments. Paired-end sequencing was performed on an Illumina NovaSeq instrument (Illumina, Inc., United States) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China) using NovaSeq reagent kits (Illumina, Inc., United States).

Metagenomic data were assembled using MEGAHIT (Li et al., 2015)¹. The predicted open reading frames with lengths ≥100 bp were retrieved and translated into amino acid sequences via the National Center for Biotechnology Information translation table. High-quality reads were aligned to non-redundant gene catalogs to calculate gene abundance with 95% identity via SOAPaligner (Li et al., 2008)². Representative sequences of the non-redundant gene catalog were aligned to the non-redundant database via Diamond³ (Buchfink et al., 2015) for taxonomic annotation. Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation was conducted against the KEGG database via Diamond.

The 47 BA standards were weighed and prepared as working standard solutions. These standards and the BA isotope standard solutions were used for analysis. The standard solution was added to 100 μL of the sample, which was mixed for 30 s and sonicated at low temperature for 30 min before being stored at −20°C for 30 min. The mixture was subsequently centrifuged at 4°C for 15 min at 13000 rcf, the supernatant was removed, and the solvent was removed by a stream of nitrogen. Then, 100 μL of 50% acetonitrile was added, and the mixture was vortexed for 30 s, sonicated at low temperature for 10 min (5°C, 40 kHz), and centrifuged at 4°C for 15 min at 13,000 rcf, after which the supernatant was removed for liquid chromatography–tandem mass spectrometry (LC–MS/MS) analysis (de Groot et al., 2020; Jung et al., 2021).

LC–MS/MS analysis of the samples was conducted on an ExionLC AD system coupled with a QTRAP^® 6,500+ mass spectrometer (Sciex, United States) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China). In brief, the samples were separated via a Waters BEH C18 column (150 × 2.1 mm, 1.7 μm) at 50°C, and the mass spectrometric data were collected with a SCIEX QTRAP 6500+ mass spectrometer equipped with an electrospray ionization (ESI) source operating in negative ion mode. The raw LC–MS/MS data were imported into SCIEX software OS, and the BA concentrations were calculated according to the standard curve.

SPSS 23.0 (New York, United States) and GraphPad Prism 5.0 (San Diego, United States) were used for statistical analyses and figure generation (Zhao and Fang, 2024). Normality and homogeneity of variance tests were conducted on the data. The experimental results are expressed as the means ± standard deviations. Multiple group comparisons were conducted via one-way analysis of variance (ANOVA) followed by a least significant difference test, Dunnett’s T3 test, or the Kruskal–Wallis test. Principal coordinate analysis (PCoA) was used to display the changes in classification with R software (v3.1.1), mixOmics, and the ade4 package. A Venn diagram was constructed with R software (v3.1.1). The KEGG pathways were predicted via PICRUSt2 v2.2.0-b and R (v3.4.10). ImageGP was used for intergroup pathway analysis. Functional difference analysis was performed via the Wilcoxon test. Receiver operating characteristic (ROC) curves were constructed, the areas under the curve (AUCs) were calculated, decision curve analysis was performed, and correlation network diagrams were generated with R software (v4.2.1). p-values <0.05 were considered to indicate statistical significance.

Compared to the Con group, the Mod group presented significant increases in body weight and liver weight (p < 0.01). Compared to those in the Mod group, the body weights and liver weights in the OCA and XFD groups were significantly lower (p < 0.01). In addition, among the four groups, there were no significant differences in the FBG or FIns levels or HOMA–IR scores (p > 0.05), as shown in Table 1.

Compared to those in the Con group, the serum and liver TC and TG levels in the Mod group were significantly greater (p < 0.01), whereas compared to those in the Mod group, the lipid levels in the OCA and XFD groups were significantly lower (p < 0.05, p < 0.01) (Figure 1A). Compared to those in the Con group, the serum ALT levels in the Mod group were significantly greater (p < 0.05), but compared to those in the Mod group, the serum ALT levels in the XFD group were significantly lower (p < 0.05). However, there was no difference in the levels of serum AST among the four groups (Figure 1B). Compared to those in the Con group, the BSH levels in the Mod group were significantly lower (p < 0.01), and compared to those in the Mod group, the BSH levels in the OCA and XFD groups were significantly greater (p < 0.05) (Figure 1C).

Compared to those in the Con group, a greater degree of hepatocyte steatosis and slight lymphocyte infiltration was observed in the livers of the Mod group. However, hepatocyte steatosis was sporadic in the OCA and XFD groups compared to that in the Mod group (Figure 2A). Moreover, the volume of adipocytes in white adipose tissue in the Mod group was greater than that in the Con group, and the OCA and XFD groups had fewer adipocytes than the Mod group (Figure 2B).

In addition, the oil red O-positive areas of adipocytes in the liver and white adipose tissue were significantly greater in the Mod group than those in the Con group (p < 0.01), and compared to those of the Mod group, the oil red O-positive areas of the OCA and XFD groups were significantly lower (p < 0.05, p < 0.01) (Figures 2C–E).

A total of 46,372,471 clean reads were obtained for a total of 18,105 species. Notably, 11,964 common species were found among the four groups. A total of 401 unique species were identified in the Con group, 616 in the Mod group, 508 in the XFD group, and 298 in the OCA group (Figure 3A). PCoA revealed differences in the beta diversity of the gut microbiota among the four groups. The Mod group was separated from the Con group, and the OCA and XFD groups were separated from the Mod group (Figure 3B).

The top 20 phyla, 30 genera, and 50 species identified in the gut microbiota are listed in Figures 3C–E. The main microorganisms were not significantly different between the OCA and XFD groups and the Mod group at the phylum and genus levels (p > 0.05). Compared to those of the top 20 species in the Con group, the abundances of Bacteroidales_bacterium, Prevotella_sp., and bacterium_0.1xD8-71 were significantly lower in the Mod group (p < 0.05 or p < 0.01), whereas the abundance of unclassified_g_Turicibacter was significantly greater (p < 0.01). Compared to that in the Mod group, the abundance of Bacteroidales_bacterium was significantly greater in the XFD group (p < 0.05), and the abundance of unclassified_g_Turicibacter was significantly lower (p < 0.05) (Figure 3F). On the basis of a linear discriminant analysis (LDA) score > 3.0, 22 species in the four groups were defined as dominant components of the gut microbiota. Notably, unclassified_g_Turicibacter, unclassified_g_Staphylococcus, Turicibacter_sp., Staphylococcus_nepalensis, and Turicibacter_sp_TS3 were significantly enriched in the Mod group (p < 0.01); Erysipelotrichaceae_bacterium, Allobaculum_mucolyticum, Clostridium_sp_CAG_557, Bifidobacterium_pseudolongum, and Eubacterium_coprostanoligenes were significantly enriched in the OCA group (p < 0.05 or p < 0.01); and Ligilactobacillus_animalis, Lactobacillus_intestinalis, Ligilactobacillus_apodemi, and Bacteroidales_bacterium were significantly enriched in the XFD group (p < 0.05), as shown in Figure 3G. According to LDA and differential analysis, Bacteroidales_bacterium may be a crucial species in the XFD group.

A comparison of the predicted gene sequences with those in the KEGG database revealed 12,615 KEGG functional annotations. The KEGG pathway enrichment analysis was performed at the third level, with an abundance of ≥0.0045. The metabolomic analysis revealed that NAFLD was related to the biosynthesis of amino acids; biosynthesis of cofactors; pyrimidine metabolism; glycolysis/gluconeogenesis; pyruvate metabolism; alanine, aspartate, and glutamate metabolism; peptidoglycan biosynthesis; glyoxylate and dicarboxylate metabolism; pantothenate and CoA biosynthesis; starch and sucrose metabolism; cysteine and methionine metabolism; butanoate metabolism; terpenoid backbone biosynthesis; arginine biosynthesis; fatty acid metabolism; glycerolipid metabolism; and thiamine metabolism (Figure 4A).

The XFD treatment affects a series of pathways that attenuate NAFLD in rats; 21 of these pathways that displayed differences in protein expression are shown in Figure 4B. These pathways include sesquiterpenoid and triterpenoid biosynthesis, regulation of lipolysis in adipocytes, vitamin B6 metabolism, riboflavin metabolism, and terpenoid backbone biosynthesis. The KEGG functions of the gut microbiota with an LDA score > 2.0 were defined as dominant, which led to the selection of 12 signaling pathways across the four groups. In the Mod group, pathways related to neurodegeneration and multiple diseases were significantly enriched in differentially expressed proteins (p < 0.05), and in the XFD group, pathways related to cysteine and methionine metabolism, lysine biosynthesis, sulfur metabolism, and monobactam biosynthesis were significantly enriched in differentially expressed proteins (p < 0.05) (Figure 4C). Moreover, 43 genes related to NAFLD were identified via KEGG analysis. These genes are related to the AMPK signaling pathway, the AMPK/p38 signaling pathway, the cdc42 signaling pathway, and the IRE1α pathway (Figure 4D).

Seventeen differentially abundant BAs were detected in the fecal samples (Figure 5A), whereas 24 differentially abundant BAs were detected in the serum and liver samples (Figures 5B,C). On the basis of the above findings, XFD normalized the levels of 16, 23, and 14 BAs in the feces, liver, and serum, respectively. Among 16, 23, and 14 BAs, XFD can normalize the levels of six common differentially abundant BAs: glycochenodeoxycholic acid, deoxycholic acid, murideoxycholic acid, lithocholic acid, 23-nordeoxycholic acid, and 3β-ursodeoxycholic acid.

Random forest analysis was conducted on the basis of the difference indices of 7 clinical indicators, 18 BAs, BSH, and Bacteroidales_bacterium to identify the 5 top biomarkers for the NAFLD diagnostic model. These indicators included 3β-ursodeoxycholic acid in feces, glycochenodeoxycholic acid in the serum, liver TGs, murideoxycholic acid in feces, and glycochenodeoxycholic acid in the liver (Figure 6A). The AUC of glycochenodeoxycholic acid in the serum and liver was 1.0 for NAFLD patients (Figure 6B). Importantly, there are more advantages than disadvantages in the use of glycochenodeoxycholic acid levels in the serum and liver as clinical biomarkers of NAFLD according to the decision curve analysis (Figure 6C). The correlation coefficient between the serum and liver levels of glycochenodeoxycholic acid was 0.799. Furthermore, there were negative correlations between BSH abundance and the levels of glycochenodeoxycholic acid in the serum and liver and between BSH abundance and the level of murideoxycholic acid in the feces (Figure 6D).