Male C57 mice aged 8 weeks and weighing between 18 and 22 g were randomly allocated into two groups: 9 mice in the high-fat group and 8 in the control group. The high-fat group was fed through 85% basal feed +15% lard +1.5% cholesterol; the low-fat group was fed through normal feed. The feed composition includes 97 g/kg moisture, 194.8 g/kg crude protein, 46 g/kg crude fat, 20 g/kg crude fiber, 56 g/kg crude ash, 13.4 g/kg calcium, 8.2 g/kg total phosphorus, 0.146 mg/kg aflatoxin B1, and < 10 CFU/g total colony count. All mice had ad libitum access to food and water throughout the study. The high-fat group received a diet high in -fat content, whereas the control group was provided with a standard normal diet. After a feeding period of 13 weeks, blood samples were collected from the eyeball, fecal samples were collected and liver tissue was obtained through laparotomy. This study received approval from the ethics committee of Kunming Medical University.

The blood collected from the eyeballs in each group at a range of 2–8°C by centrifugation at 3000 rpm for 15 min. Subsequently, the resulting and the supernatant underwent analysis using a biochemical analyzer (SMT-120VP, Smart).

The liver tissues obtained from the mice in each group were initially fixed in paraformaldehyde, and subsequently treated with the OCT embedding agent, and rapidly frozen and embedded on the frozen table of microtome. To visualize lipid content, frozen sections were immersed in an oil red O working solution in darkness for 8–10 min. Afterwards, differentiation was achieved by briefly submerging the sections in 60% isopropyl alcohol for 3–5 s. The staining of nuclei was accomplished using a haematoxylinstaining solution (G1004, Servicebio), followed by a series of washing steps with water to remove excess stain. Once washed, dried slightly, the sections were sealed using a glycerine gelatine sealing medium to preserve the samples.

The total fecal microbiota DNA was extracted and the V3-V4 regions of the 16S rRNA gene were amplified. DNA extraction kits were available for different types of samples in order to ensure DNA extraction efficiency and quality. (OMEGASoilDNAKit; OMEGAWaterDNAKi; OMEGAStool DNAKit). The PCR primer (Supplementary Table 1) was designed against the conserved region to target the variable region of the 16S /ITS2 rDNA gene. After 35 cycles of PCR, sequencing adapters and barcodes were added for amplification. PCR amplification products were detected by 1.5% agarose gel electrophoresis. The target fragments were recovered using the AxyPrep PCR Cleanup Kit. The PCR product was further purified using the Quant-iT PicoGreen dsDNA Assay Kit. The library was quantified on the Promega QuantiFluor fluorescence quantification system. The pooled library was loaded on Illumina platform using a paired-end sequencing protocol (2 × 250 bp). The composition of the fecal microbiota was evaluated using the quantitative insights into microbial ecology (QIIME) software (Release 138) for microbiota analysis.

Additionally, the collected tissue samples were thawed on ice, and metabolites were extracted using a 50% methanol buffer. Briefly, 20 μL of each sample was mixed with 120 μL of precooled 50% methanol, vortexed for 1 min, and then incubated at room temperature for 10 min. The extraction mixture was stored overnight at −20°C. The next day, the mixture was centrifuged at 4,000 g for 20 min, and the supernatants were carefully transferred into new 96-well plates. The samples were then stored at −80°C until further analysis by LC–MS. In addition, quality control (QC) samples were prepared by combining 10 μL of each extraction mixture. QC samples with a relative standard deviation (RSD) greater than 30% were removed from the analysis to ensure data quality. Subsequently, the extracted metabolites were annotated and quantified using the metaX software in combination with the KEGG and HMDB databases to identify and characterize the metabolites of interest.

The sequencing depth was assessed using the Good’s coverage index. Various a Alpha-diversity indexes, including Chao1, observed species, Shannon and Simpson, were compared to analyze the species richness and uniformity (Mohakud et al., 2019). Beta-diversity analysis was used to evaluate the species’ complexity among different groups. In this study, principal coordinate analysis was used to illustrate the differences between the different groups. Subsequently, the analysis of similarities (Anosim) was used to analyze whether differences between the two groups were significantly greater than those within the groups (Draiko et al., 2019).

The GM composition was further analyzed, comparing the distribution of the GM at both the phylum and genus levels. Characteristic strains were identified using linear discriminant analysis effect size (LEfSE), and operational taxonomic units with a relative abundance exceeding 0.5% were selected for analysis based on the linear discriminant analysis (LDA) scores (Segata et al., 2011).

Receiver operating characteristic (ROC) curves were calculated using the “pROC” R package (version 1.18.0), and the area under the ROC curve (AUC) values were used to access the distinguish ability of characteristic strains (Robin et al., 2011). Functional prediction of the diagnostic strains (AUC > 0.8) was assessed, and differences between NAFLD mice and HCs were studied using the “PICRUSt2” R package (version 0.2.3) (Douglas et al., 2020).

The metabolic differences between the NAFLD mice and HCs were compared through the construction of an orthogonal partial least squares discriminant analysis (OPLS-DA) model. This model’s validity was assessed using a permutation test. Subsequently, three methods, were employed to select metabolites: the multiple change method (FC value), the t test method (p value, and FDR value), and the partial least squares discriminant analysis method (VIP value). The differentially expressed metabolites between the two groups were identified using the following criteria: VIP > 1, FC ≥ 1, and p < 0.05. In addition, enrichment analysis of these differentially expressed metabolites was performed using the online platform “Metaboanalyst” (Chong et al., 2019).

The mRNA expression levels between the NAFLD mice and HCs were compared using the “limma” R package (version 3.52.4) (|logFC| > 0.5, p < 0.05), and visualized using the ‘ggplot2’ R package and Kolde R (2019). _pheatmap: Pretty Heatmaps_. R package version 1.0.12.¹ Then, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEGs were conducted using the ‘clusterprofiler’ R package (version 4.4.4) (Wu et al., 2021).

The correlations between the characteristic strains and differential metabolites were assessed using the ‘Spearman’ package (|r| ≥ 0.3, p < 0.05). Metabolites in the pair were regarded as the key metabolites 1, and diagnostic strains were regarded as key microbiota. Similarly, the correlations between the DEGs and differential metabolites were studied, and the DEGs in the pair were regarded as key genes, whereas metabolites were regarded as key metabolites 2. Collinear networks of the key microbiota-metabolites 1 and key genes-metabolites 2 were constructed using ‘Cytoscape’ (version 3.7.1) (Shannon et al., 2003). Finally, key metabolites 1 and 2 were combined, and a collinear network of key genes-microbiota-metabolites was constructed using ‘Cytoscape.’

Total RNA was extracted from the control and model groups, and the reverse transcription reaction was subsequently performed, followed by qPCR reaction with 2 × Universal Blue SYBR Green qPCR Master Mix kit. The PCR cycling conditions were as follows: pre-denaturation at 95°C for 1 min; denaturation for 20 s at 95°C, annealing for 20 s at 55°C,and final extension for 30 s at 72°C; total, 40 cycles.

Statistical analysis was performed using R software (version 4.2.3) and the maps were plotted using Prism GraphPad, and the experimental data are presented as mean ± SEM. The Welch’s t-test was used to compare the measurement data between groups. p < 0.05 for statistical significance; p < 0.01 for high statistical significance; p < 0.001 for extremely high statistical significance.

Upon examination, liver tissues from the model group of C57 mice exhibited a distinctive yellow hue and enlarged sizes (Figure 1A). Moreover, the levels of TG, TC and LDL were significantly increased, whereas HDL levels showed a significant decrease (Figure 1B). These findings collectively suggest the successful establishment of the NAFLD model (Figure 1).

Excessive fat deposition represents a crucial pathological characteristic of NAFLD. To validate the successful establishment of the model, we conducted observations of adipogenesis using oil red O staining. Results demonstrated a substantial increase in fat content within the model group compared to the control group, further substantiating the successful construction of the model (Figure 2).

The goods coverage curve approaching 1, indicated that sufficient sequencing data were acquired, affirming a reasonable sequencing depth (Figure 3A). No significant differences in alpha-diversity indexes were observed between the NAFLD and HC groups, suggesting that the richness and uniformity of bacterial communities remained largely unchanged owing to NAFLD (Figure 3B). However, the PCoA results showed a distinct clustering of patients with NAFLD and HCs (Figure 3C). Furthermore, the ANOSIM analysis demonstrated greater differences between groups than within groups (ANOSIM statistics = 0.6233, p = 0.001), signifying the meaningfulness of grouping and substantial alterations in the structure and composition of the GM induced by NAFLD.

At the phylum level, the gut microbiota primarily comprised Firmicutes, Bacteroidetes, Actinobacteria, Desulfobacterota, Patescibacteria, Proteobacteria, Verrucomicrobiota, Campylobacterota, Deferribacterota, and Synergistota. Comparing the NAFLD and HC groups, revealed significant decreases in the abundances of Firmicutes and Synergistota in the NAFLD groups (p < 0.05), whereas Patescibacteria and Deferribacterota showed decreased abundance, Others demonstrated significantly decreased abundances in the NAFLD group (p < 0.05) (Figure 4A and Supplementary Table 2). At the genus level, the NAFLD group exhibited prominent representation of Muribaculaceae, Allobaculum, Dubosiella, Ligilactobacillus, Ruminococcaceae, Paramuribaculum, Bifidobacterium, and Desulfovibrio, among others. The HC group showed a predominant presence of Ligilactobacillus, Muribaculaceae, Lactobacillus, HT002, Dubosiella, Allobaculum, Candidatus, and Clostridia, among others. Notably, the NAFLD group demonstrated significantly increased abundances of Muribaculaceae, Allobaculum, Paramuribaculum, and Ruminococcaceae, and decreased abundances of Ligilactobacillus, Lactobacillus, and HT002 (Figure 4B and Supplementary Table 3). Further investigation using LEfSE analysis identified 46 differentially abundant key gut microbiota, including 35 with notable roles in the NAFLD groups and 11 in the HC groups. Among these, Ileibacterium, Ruminococcaceae (Firmicutes), Olsenella, Duncaniella and Paramuribaculum (Bacteroidota), and Bifidobacterium, Coriobacteriaceae_UCG_002 and Olsenella (Actinobacteriota) were characteristic strains of NAFLD. On the other hand, HT002, Lactobacillus (Firmicutes), Bacteroides (Bacteroidota) were characteristic strains of HC (Figures 4C,D and Supplementary Table 4).

In this study, the ROC curves of the characteristic strains were calculated, and 13 diagnostic strains were obtained with an AUC value >0.8 (AUC_{Bifidobacterium} = 0.969, AUC_Olsenella = 1.000, AUC_{Enterorhabdus} = 0.812, AUC_Bacteroides = 0.812, AUC_Duncaniella = 1.000, AUC_Muribaculum = 0.922, AUC_{Paramuribaculum} = 0.953, AUC_{Ileibacterium} = 0.891, AUC_Turicibacter = 0.922, AUC_HT002 = 0.883, AUC_{Lactobacillus} = 0.979, AUC_{Ruminococcaceae} = 0.984, AUC_{Ileibacterium} = 0.891) (Figure 5A). In addition, the functional prediction of the diagnostic strains was assessed, and the results showed that 121 pathways were significantly different between the two groups. Among these groups, methylerythritol phosphate pathway, glycogen degradation, l-ornithine, l-valine, l-tryptophan, and l-isoleucine biosynthesis pathways, among others, were significantly highly expressed in the NAFLD groups (p < 0.05) (Figure 5B).

The constructed OPLS-DA model revealed distinct differences in metabolite composition between the two groups (Q2 = 0.837, R2Y = 0.999) (Figure 6A). A total of 2,497 differentially regulated metabolites, including 947 down-regulated and 1,551 up-regulated metabolites, were identified in the NAFLD group, meeting the criteria of VIP > 1, FC ≥ 1, and p < 0.05 (Figures 6B,C). In addition, these differentially regulated metabolites were predominantly enriched in pathways related to tryptophan, alpha-linolenic acid, linoleic acid and methylhistidine metabolism (Figure 6D).

A total of 2,510 DEGs (1,531 down-regulated and 979 up-regulated) were screened in the NAFLD groups compared with the HC groups (Figures 7A,B). In addition, these DEGs exhibited significant enrichment in 92 GO, categories, notably involving processes such as small molecule catabolic, processes and the regulation of small GTPase mediated signal transduction (Figure 7C). Furthermore, the KEGG results revealed associations of these DEGs with pathways related to cholesterol metabolism, amino acid biosynthesis, and primary bile acid biosynthesis, among others (Figure 7D).

As shown in Figure 8A, characteristic strains of the HC groups—Bacteroides, HT002 and Lactobacillus demonstrate strong negative correlations with docosahexaenoic acid, n-isobutyrylglycine and taurine. Conversely, the key microbiota of the NAFLD groups exhibit strong positive correlations with cholanoic acid, 5-methyltetrahydrofolic acid etc. The correlation between the key metabolites 2 and genes is shown in Figure 8B. Shisa8 and Lrrc8e display strong positive correlations with cis-4-hydroxyequol, whereas Gm38220 and Gm42977 exhibit strong negative correlations with Glu-Ile. Moreover, the constructed network involving 8 key strains, 14 key genes and 83 key metabolites elucidates potential regulatory interactions. Notably, Mup-ps16 has been indicated to potentially regulate the metabolism of 3,7-dihydroxy-12-oxocholanoic acid through interactions with Duncaniella and Olsenella. Additionally, Proser2 shows potential regulatory capabilities across the metabolism influenced by most diagnostic microbiota (Figure 8C and Supplementary Table 5).

Based on the results of multi-omics joint analyses, we further verified the expression of differential genes using RT-PCR. The results showed that compared with that in the control group, the expression of Slc22a7 in the model group was significantly increased, whereas the expressions of Hsd3b5, Zfp334, Ace2, Dbp and Proser2 were significantly decreased. Cyp2b9 and Ccr7 showed no significant differences. In addition to Dbp, Hsd3b5, Cyp2b9 and Ccr7, Slc22a7, Zfp334, Ace2 and Proser2 expressions were consistent with the results of multi-omics analysis (Figure 9).