Lingguizhugan decoction (powder batch number: Z201101) was provided by Professor Tong Zhang, School of Pharmacy, Shanghai University of Traditional Chinese Medicine. Poriacocos (Schw.) Wolf (batch number: Y2003002), Cinnamomum cassia Presl (batch number: 200608), Atractylodes macrocephala Koidz. (batch number: YP200601), and Glycyrrhiza uralensis Fisch. (batch number: YP200601) in a 2:1.5:1:1 ratio, added 12 times the amount of water, decoction 2 times, every 1.5 hours, collected the first water decoction, and set aside. Decoction filtered, combined two filtrates, the filtrate concentrated to a relative density of 1.07 ~ 1.09 (65 ± 5°C), spray-dried, and crushed into a fine powder for use (4.56 g crude medicine extracted 1 g extract powder, which was stored in a dry environment). All herbs were purchased from Jiangsu Sanhexing Chinese Medicine Research Co., Ltd. (Jiangsu, China). The fingerprint was used to control the quality of the LGZG ( Supplementary Figure 1 ).

Male Sprague-Dawley rats (8-week-old, 300-350 g), were purchased from Shanghai Jihui Experimental Animal Technology Co., Ltd. (Shanghai, China), and maintained in the Laboratory Animal Center of Shanghai University of Traditional Chinese Medicine, which were kept in a constant temperature and humidity with a 12 h light/dark cycle. After one week of adaptive feeding, the rats were randomly divided into two groups. The Control group (n = 6) was fed a chow diet, and the model group was fed an MCD diet (Lot No.: 21011101, purchased from Changzhou SYSE Bio-Tech. Co., Ltd. (Changzhou, China). After two weeks, the model group rats were randomly divided into two groups: the NASH group and the Lingguizhugan decoction-intervened (LGZG) group (n = 6 for each group). The NASH group was fed the MCD diet, and the LGZG group was fed the MCD diet supplemented with LGZG, which was administered at 0.1 mL/10g body weight (clinical equivalent dose:16.56 g crude drug/kg) by oral gavage every day for four weeks. The Control group and NASH group rats were given an equal volume of pure sterilized water. Body weight and food intake were recorded every two days to observe the state of the rats. At the end of the experiment, the rats were fasted for 12 h, weighed, and injected with 10% chloral hydrate for anesthesia. Blood was then taken from the abdominal aorta and the serum was separated for biochemical analyses, while the liver was removed, weighed, and repacked. In addition, the ileum, colon, and feces of the ileocecal region were collected. All animal experiment procedures were approved by Shanghai University of Traditional Chinese Medicine Animal Experiment Ethics Committee (No. PZSHUTCM201218010).

Serum alanine aminotransferase (ALT, batch number: 01ALT210107), and aspartate aminotransferase (AST, batch number: 01AST201111) levels were analyzed and detected by an automatic biochemical analyzer (TBA-40FR, Toshiba, Tokyo, Japan). The supernatant of liver homogenate was prepared by freezing liver tissue. Then the levels of liver triglyceride (TG, batch number: 02TG201209), and liver total cholesterol (TC, batch number: 01CHOL201030) were measured using an automatic biochemical analyzer. All kits were purchased from Shanghai Huachen Biological Reagent Co., Ltd. (Shanghai, China).

The left upper half of the liver tissues of rats were fixed in 10% neutral formalin for one week, then dehydrated, paraffin-embedded, and sectioned into 4 μm slices. Pathological changes of the liver tissues were evaluated by hematoxylin and eosin (H&E, BaSO, China) staining and the total NAFLD Activity Score (NAS score), as previously reported (Kleiner et al., 2005). The scoring system was mainly comprised of three histological features, which were evaluated semi-quantitatively: steatosis (0-3), lobular inflammation (0-3), and hepatocellular ballooning (0-2). Liver tissues used for Oil Red O (ORO) staining were frozen and embedded in optimal cutting temperature compound (OCT, SAKURA Tissue-Te, America), which were sectioned into 10 µm slices, stained with ORO and counter-stained with hematoxylin. The areas of liver ORO staining were morphologically analyzed using a StrataFAXS II image analysis system (Strata FAXS II, Vienna, Austria).

Hepatic levels of proinflammatory cytokines such as tumor necrosis factor α (TNF-α, batch number: ml002859-2), interleukin 6 (IL-6, batch number: ml102828-2), and interleukin 1β (IL-1β, batch number: ml003057-2) were measured using enzyme-linked immunosorbent assay (ELISA) kit (Shanghai Enzyme-linked Biotechnology Co., Ltd, Shanghai, China), according to the manufacturer’s protocol.

The 16S rRNA of fecal samples was sequenced by Metabo-Profile Biotechnology (Shanghai) Co., Ltd. (Shanghai, China), following the manufacturer’s procedures. Briefly, total genomic DNA was extracted using the OMEGA Soil DNA Kit (M5635-02) (Omega Bio-Tek, Norcross, GA, USA). The quantity and quality of extracted DNAs were measured using a NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis, respectively. PCR amplification of the V3–V4 regions of bacterial 16S rRNA genes was then conducted. PCR amplicons were purified with Vazyme VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China) and quantified using the Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). After the individual quantification step, amplicons were pooled in equal amounts, and paired-end 2×250 bp sequencing was performed using the Illumina NovaSeq platform with the NovaSeq 6000 SP Reagent Kit (500 cycles; Illumina, San Diego, CA, USA).

The sequence data were processed using QIIME2 according to previously described methods (Gupta et al., 2022). Briefly, raw sequence data were demultiplexed, quality filtered, denoised, merged, and chimera removed. Non-singleton amplicon sequence variants (ASVs) were aligned to construct a phylogeny. Taxonomy was assigned to ASVs using the classify-sklearn naiïve Bayes taxonomy classifier against the SILVA Release 138 Database (Gupta et al., 2022). Alpha diversity indices, including the Chao1 richness estimator, Shannon diversity index, and Simpson index were calculated and visualized as box plots. Un-weighed principal coordinate analysis (PCoA) and the unweighted pair-group method with arithmetic means (UPGMA) system clustering tree for beta diversity were performed to investigate the structural variation and similarity of microbial communities across samples. The relative abundance of taxonomy at the genus level was visualized as a bar chart using R software (R Foundation for Statistical Computing, Vienna, Austria). Metastats analysis with p values less than 0.05 was used to identify significantly changed microbiota between groups. A Venn diagram was performed to obtain overlapped microbiota between pairwise compared groups (NASH vs. Control and LGZG vs. NASH), and a hierarchical cluster was used to reveal changing patterns of microbiota among groups and to identify LGZG reversed microbiota.

Metabolomics analyses of fecal samples were performed by Metabo-Profile Biotechnology (Shanghai) Co., Ltd. (Shanghai, China), following the manufacturer’s instructions. Feces samples were thawed on an ice bath, about 5mg of each lyophilized sample was weighed and transferred to a new 1.5 mL tube. Then, 25 μL of water was added, and the sample was homogenized with zirconium oxide beads for 3 minutes, then 120 μL of methanol containing an internal standard was added to extract the metabolites. The sample was homogenized for another 3 minutes and then centrifuged at 18000 g for 20 minutes. A total of 20 μL of supernatant was transferred to a 96-well plate. The following procedures were performed using an Eppendorf epMotion Workstation (Eppendorf Inc., Hamburg, Germany). A total of 20 μL of freshly prepared derivative reagents was then added to each well. The plate was sealed and derivatization was conducted at 30°C for 60 min. After derivatization, 330 μL of an ice-cold 50% methanol solution was added to dilute the sample. Then the plate was stored at -20°C for 20 minutes, followed by 4000g centrifugation at 4°C for 30 minutes. A total of 135 μL of supernatant was then transferred to a new 96-well plate with 10 μL internal standards in each well. Serial dilutions of derivatized stock standard were added to the left wells.

Ultra-performance liquid chromatography coupled to tandem mass spectrometry (UPLC-MS/MS) system (ACQUITY UPLC-Xevo TQ-S, Waters Corp., Milford, MA, USA) was used to quantify all targeted metabolites. The following were the instrument settings. Column: ACQUITY UPLC BEH C18 1.7 µM Van Guard pre-column (2.1×5 mm) and ACQUITY UPLC BEH C18 1.7 µM analytical column (2.1×100 mm); column temperature is 40°C. The mobile phase consisted of distilled water containing 0.1% formic acid (A) and 70% acetonitrile-30% isopropanol (B). The following gradient conditions were used: 0-1.0 min, 5% B; 1.0-11.0 min, 5% B - 78% B; 11.0-13.5 min, 78% B - 95% B; 13.5-14.0 min, 95% B - 100% B; 14.0-16.0 min, 100% B; 16.0-16.1 min, 100% B - 5% B; 16.1-18.0 min, 5% B; flow rate (mL/min):0.40; and injection volume (µL):5. Additionally, the analytical quality control experiments were performed according to previously described methods (Xie et al., 2021). The pooled QC samples were prepared by mixing aliquots of the study samples such that the pooled samples broadly represented the biological average of the whole sample set. The QC samples for this project were prepared with the test samples and injected at regular intervals (after every 10-14 test samples for LC-MS) throughout the analytical run to ensure a consistently high quality of analytical results. All the QC samples, calibrators, and blank samples were analyzed across the entire sample set to diminish analytical bias.

The raw data files generated by UPLC-MS were processed using the MassLynx software (v4.1, Waters, Milford, MA, USA) to perform peak integration, calibration, and quantitation for each metabolite. Multivariate statistical analyses, including principal component analysis (PCA) and orthogonal partial least square discriminant analysis (OPLS-DA), were performed using SIMCA software (version 14.1). Besides, a 200-times permutation was conducted to validate the OPLS-DA model against over-fitting, and variable influence on projection (VIP) values of each metabolite in the OPLS-DA model was calculated. For univariate statistical analyses, a t-test or Wilcoxon test was performed to calculate p values. Finally, VIP>1 and p<0.05 were used to identify significantly changed metabolites between groups. A Venn diagram was performed to obtain overlapped metabolites between pairwise compared groups (NASH vs. Control and LGZG vs. NASH), and hierarchical clustering was used to identify the changing patterns of metabolites among groups and to identify LGZG reversed metabolites. Enrichment analyses of reversed metabolites were performed using MetaboAnalyst online tools.

The measurement data were shown as the mean ± standard deviation (SD), and all data were statistically analyzed by one-way analysis of variance (ANOVA) followed by independent-sample-t-test using GraphPad Prism 8.4.2 (GraphPad, San Diego, CA, USA). Spearman rank correlation was performed among LGZG reversed gut microbiota, fecal metabolites, and NASH indexes. A value of p <0.05 was considered statistically significant.

To evaluate the effect of LGZG on NASH, we constructed the MCD diet-induced NASH model. MCD diet-fed rats were then treated with either LGZG or vehicle for four weeks. Compared to the Control group, the NASH group showed a significant decrease in body weight and an increase in liver index (the ratio of liver weight to body weight), serum ALT and AST levels. In addition, MCD feeding significantly elevated liver TG content while decreasing liver TC content in rats. After four weeks of LGZG intervention, serum AST, ALT, and liver TG levels were significantly decreased in MCD diet-fed rats. However, body weight and liver index were not significantly changed after the decoction intervention ( Figure 1 ).

To observe the effects of LGZG on liver histological changes in NASH, a histopathological examination was performed. As shown in Figure 2 , H&E and ORO staining revealed that steatosis, inflammation, and lipid droplets were increased in the MCD diet-induced NASH rats compared with the Control group, which were partially restored by LGZG. Moreover, the inflammation score and NAS score were significantly higher in the NASH than in the Control group, which were significantly lowered by LGZG intervention. The ORO staining areas of the LGZG group were reduced by nearly 1/3 compared with the NASH group (average 64.8% to 43.1%). Meanwhile, the level of liver proinflammatory cytokine IL-1β was significantly increased in NASH rats, while IL-1β, IL-6, and TNF-α levels were significantly decreased by the LGZG.

The 16S rRNA sequencing was used to investigate the effect of LGZG on gut microbiota in MCD diet-induced NASH rats. The α-diversity indexes including Chao1, Shannon, and Simpson were used to determine the ecological diversity within the microbial community. As shown in Figure 3A , there were decreased trends of Chao1, Shannon, and Simpson in the NASH group compared with the Control group; however, there was no significant difference among the Control, NASH, and LGZG groups. To investigate structural variations of microbial communities across samples, PCoA of β-diversity was performed. The results revealed that there were distinct separations, which indicated that there might be different compositions of gut microbiota among the control, NASH, and LGZG groups ( Figure 3B ). Furthermore, UPGMA analysis revealed that the composition of gut microbiota was similar within the same group and different among groups ( Figure 3C ). To reveal the different compositions, the relative abundances of gut microbiota at the genus level were compared across samples in the three groups. As shown in Figure 3D , the relative abundances of many microbiotas at the genus level were different among the Control, NASH, and LGZG groups. For example, the relative abundance of Lactobacillus accounted for more than 40% of gut microbiota in the Control group, which decreased to about 5% in the NASH group and was restored to about 10% in the LGZG group.

To identify the significantly changed microbiota at the genus level among groups, Metastats analyses with values of p < 0.05 were performed. The results revealed 26 significantly changed genera (such as Ruminococcus_1, Peptococcus, Lactococcus, etc.) between the NASH and Control groups, and 51 significantly changed genera (such as Eisenbergiella, Intestinimonas, Allobaculum, etc.) between the LGZG and NASH groups. The 10 representatives significantly changed genera were shown in Figures 4A, B , and the detailed information of all significantly changed genera was listed in Tables 1 , 2 , and Supplemental Table 1 . Besides, a Venn diagram was used to identify overlapped genera between the pairwise groups (NASH vs. Control and LGZG vs. NASH), with a total of 13 genera being identified ( Figure 4C ). Furthermore, hierarchical clustering was then used to identify the microbial changed patterns among groups and to identify LGZG reversed genera. As shown in Figure 4D , 11 genera exhibited an opposite pattern between NASH vs. Control and LGZG vs. NASH, including Ruminococcus_1, A2, Odoribacter, Harryflintia, Peptococcus, Ruminococcaceae_UCG-004, Butyricimonas, Prevotellaceae_Ga6A1_ group, Prevotellaceae_UCG-001, Rikenellaceae_RC9_gut_group, and uncultured_ bacterium_o_Rhodospirillales. For example, the relative abundances of Ruminococcus _1, Odoribacter, and Butyricimonas were increased in NASH, which was decreased after LGZG treatment. Together, the results indicated that LGZG might partially reverse the change of 11 genera in NASH.

UPLC-MS was used to profile fecal metabolites across the Control, NASH, and LGZG -intervened groups. Multivariate and univariate statistical analyses were performed to identify significantly changed metabolites among groups. As shown in Figures 5A, B , both PCA and OPLS-DA models revealed that there were distinct separations among the Control, NASH, and LGZG groups, which indicated that profiles of metabolites exhibited different patterns among the three groups. In addition, PCA and OPLS-DA models for two pairwise groups (NASH vs. Control and LGZG vs. NASH) also revealed good separations ( Figures 5C–F ). To validate the OPLS-DA model against overfitting, two hundred times permutation tests were performed. As shown in Figure 6A , the results suggested good reliability of the OPLS-DA model for the Control and NASH groups with R2= (0.0,0.75), and Q2= (0.0, -0.174). Moreover, the results ( Figure 6B ) also suggested the good quality of the OPLS-DA model for the LGZG and NASH groups with R2= (0.0,0.746), and Q2= (0.0, -0.364). Furthermore, VIP values of each metabolite were obtained in OPLS-DA models ( Figures 6C, D ).

With VIP>1 and p <0.05, we obtained a total of 99 significantly changed metabolites (such as tyrosine, pyruvic acid, TDCA, etc.) between the NASH and Control groups, and 109 significantly changed metabolites (such as valeric acid, TCA, valine, etc.) between the LGZG and NASH groups. The detailed information was listed in Supplemental Tables 2 and 3 , respectively. Using a Venn diagram ( Figure 7A ), 58 overlapped metabolites (such as DCA, TDCA, isocaproic acid, etc.) were obtained between pairwise groups (NASH vs. Control and LGZG vs. NASH), with the detailed information listed in Table 3 . Moreover, hierarchical clustering revealed that 57 metabolites exhibited opposite patterns between pairwise groups ( Figure 7B ). For example, DCA and TDCA levels in NASH rats were higher than in the control group but were lowered by LGZG- intervention. Isocaproic acid was lower in the NASH group compared with the Control group, but it was increased by LGZG decoction intervention. The results indicated that LGZG decoction might reverse the alteration of metabolites in MCD diet-induced NASH. Furthermore, the 57 reversed metabolites were classified and subjected to enrichment analysis. The results revealed that most of the metabolites were amino acids, bile acids, and organic acids ( Figure 7C , Supplemental Table 5 ), which were enriched in many Kyoto Encyclopedia and Genes and Genomes (KEGG) pathways such as glutathione metabolism, glucose-alanine cycle, glycine and serine metabolism, etc. ( Figure 7D ).

To examine the relationships among LGZG-restored gut microbiota, fecal metabolites, and NASH indexes, Spearman’s correlation analyses were performed. The results revealed that most of the LGZG restored gut microbiota and its metabolites were significantly correlated with NASH indexes ( Figures 8A, B ). For example, the fecal metabolites oleylcarnitine and isoDCA were positively correlated with liver TG, NAS score, serum AST and ALT, and TDCA was positively correlated with liver TG, NAS score, and serum ALT; whereas, valine and alanine were negatively correlated with liver TG, NAS score, serum AST and ALT. In addition, the gut microbiota Prevotellaceae_ Ga6A1_group, A2, and Butyricimonas were positively correlated with the NAS score, serum AST and ALT. Interestingly, we also observed that most of the LGZG-restored gut microbiota were significantly correlated with LGZG -restored fecal metabolites ( Figure 8C ). For example, Ruminococcus_1, Odoribacter, and Peptococcus were positively correlated with some bile acids such as DCA, TDCA, TCA, etc., but negatively correlated with other fecal metabolites such as alanine, asparagine, glutamic acid, etc. The results indicated that LGZG might improve NASH partially by modulating gut microbiota and correlated metabolites.