Characteristics of gut microbiota in patients with metabolic associated fatty liver disease

The clinical and physical variables

Table 1 presents the general characteristics of the 106 included participants. The mean age was approximately 36 years in both healthy and MAFLD groups. Of the 106 participants, 81 (68 men, 13 women) were in the MAFLD group. Body mass index (BMI) level was significantly lower in the healthy volunteer group compared with the MAFLD group. Table 2 presents the general information of the MAFLD subgroups, with or without liver enzyme abnormality. No significant difference was observed in BMI between the two subgroups (P = 0.057).

We next analyzed the differences in the physical variables between the two groups. The results shown that the levels of gamma-glutamyl transferase (GGT; Fig. 1A), alanine aminotransferase (ALT; Fig. 1B), and aspartate aminotransferase (AST; Fig. 1C) were significantly higher in subjects with MAFLD. From the analysis of MAFLD subgroups, we also found that the levels of clinical variables including controlled attenuation parameter (CAP; Fig. 1D), liver stiffness measurement (LSM; Fig. 1E), fasting insulin, hypersensitive C-reactive protein, lactate dehydrogenase (LDH), triglyceride (TG), glucose (GLU), type III procollagen peptide N terminal (PIIINP), and collagen type IV (C-IV) were significantly higher in MAFLD subjects with abnormal liver enzyme but not in subjects with normal liver enzyme (Table 2).

The clinicopathological indicators involved in patients. (A-C) The boxplot shows GGT, ALT and AST of healthy and MAFLD subjects. (D-E) Significantly different clinical parameters are presented in MAFLD subgroups with or without liver enzyme abnormal. **P < 0.01, ****P < 0.0001.

Analysis of gut microbiota using metagenomic data in species level

To gain insight into microbiome diversity, Shannon, Simpson, and InvSimpson diversity indices were determined for alpha diversity analysis of the metagenome data. The Shannon Simpson, and InvSimpson diversity indexes indicated that no significant difference of microbiome diversity was observed in the healthy volunteers and MAFLD groups (Fig. 2A–C). Similarly, no significant difference existed among MAFLD subjects with abnormal liver enzyme and with normal liver enzyme (Fig. 2D,E). But we observed a tendency of increase of the microbiome diversity in the health volunteer group and MAFLD subjects with abnormal liver enzyme.

Diversity of the gut microbiome characteristics base on metagenomic sequencing in species level. (A–C) Alpha diversity indices of the intestinal bacterial communities of healthy control and MAFLD group. (D–F) Alpha diversity indices of the intestinal bacterial communities of different MAFLD subgroups. Shannon, InvSimpson, and Simpson indices. (G,H) Different letters above the bars indicate significant differences as determined by PERMANOVA. Samples are identified by filled circles. In both PCoA1 and PCoA2 principal components (PCoA1 and PCo2) are plotted. The percentage of variance in the dataset explained by each axis is reported. (I,J) LefSe analysis showed predominant gut microbiota. *P < 0.05, ns P > 0.05.

Principal coordinate analysis provided an overview of the gut microbiome and reflected the β-diversities of the different groups. The β-diversity was clearly higher in the healthy volunteer group than that in the MAFLD group (Fig. 2G). However, the β-diversity showed no difference among the subgroups of MAFLD subjects with abnormal and normal liver enzyme (Fig. 2H). To identify the predominant gut microbiota in difference group, LEfSe analysis was performed. The results showed that there were a number of different genera of gut microbiota between the healthy and the MAFLD groups, and a trend could be observed that in the MAFLD groups the relative richness of Bacteroides vulgatus was much higher than that in the healthy group. Moreover, Ruminococcus Obeum and Alistipes putredinis were highly enriched in the healthy group (Fig. 2I). We also found that Bacteroides caccae were highly enriched in normal liver enzyme group (Fig. 2J).

Analysis of gut microbiota using 16S data in genus level

Microbial 16S rRNA gene amplicons were denoised into amplicon sequence variants (ASVs) to provide a high-resolution^19,20 view of how MAFLD impacted the community structure. To uncover the microbiota distribution and genera in MAFLD and healthy groups, fecal samples were performed 16S rRNA sequencing and Shannon, Simpson, and InvSimpson indexes calculated. No diversities of microbiome were significantly different between MAFLD and healthy group (Fig. 3A–C). Subgroups were derived according to different liver enzyme for further analysis. The result shown that the microbiome diversities did not differ significantly between subgroups (Fig. 3D–F).

Diversity of the gut microbiome characteristics base on 16S data in genus level. (A–C) Alpha diversity indices of the intestinal bacterial communities of healthy control and MAFLD group. (D–F) Alpha diversity indices of the intestinal bacterial communities of different MAFLD subgroups. Shannon, InvSimpson, and Simpson indices. (G,H) Different letters above the bars indicate significant differences as determined by PERMANOVA. Samples are identified by filled circles. In both PCoA1 and PCoA2 principal components (PCoA1 and PCoA2) are plotted. The percentage of variance in the dataset explained by each axis is reported. (I,J) LefSe analysis showed predominant gut microbiota. **P < 0.01, ns P > 0.05.

To display similarities of microbiome compositions between difference group, β-diversity was calculated using UniFrac method, and PCoA was performed. The results presented a significantly different distribution among healthy and MAFLD groups using permutational multivariate analysis of variance analysis (PERMANOVA) (Fig. 3G). But did not differ between MAFLD subgroups (Fig. 3H). These results suggest that the MAFLD group had significant different gut microbiome compositions from healthy group. Meanwhile, with the aid of the LEfSe methods, we conducted analysis of different species between groups and identified a series of bio- markers as shown in the Fig. 3I,J. The resulting cladogram presents the structure of the gut microbiome and the predominant bacteria in the healthy and MAFLD groups. Notably, results of LEfse analysis also showed that Alistipes higher enriched in MAFLD subgroup of liver enzyme normal compared with liver enzyme abnormal subgroup. Alistipes belongs to the family Rikenellaceae, which was also decreased in NAFLD patients based on a study by Zhu, L. et al.²¹. And several studies have linked the presence of Alistipes genus with a healthy state²². Therefore, the restored intestinal Alistipes communities contributed to the ameliorated MAFLD.

Microbe‐set enrichment analysis (MSEA) of gut microbiota

MSEA analysis²³ can be incorporated with microbiome profiling pipelines to determine the mechanisms underlying host‐microbiome interactions. Then, we analyzed the 16S microbiome profiling data with a focus on characterizing the functions of microbes with differential abundance (DA) in MAFLD compared to healthy controls using MSEA. We prioritized human genes enriched for those MAFLD-related DA microbes. Among the top enriched genes, we found that several interesting microbe-related-gene associations (Fig. 4A,B). For example, Dorea, Lactobacillus, Megasphaera are enriched in MAFLD group. We next performed enrichment analysis for the top genes that are enriched from the DA microbes in MAFLD compared to healthy controls. The results showed that most genes were enriched in MAFLD and Hepatitis C pathway in MAFLD group (Fig. 4C,D).

Microbe‐set enrichment analysis (MSEA) of fecal microbiota. (A-B) Bipartite graph visualizing the enriched mammalian genes with their associated gut microbes in healthy controls versus MAFLD. Genes are charted as orange round nodes whereas gut microbes are plotted as yellow square. The sizes of the nodes are proportional to the number of edges in the bipartite graph whereas the width of the edges indicates the strength of the enrichment measured by combined scores from the MSEA algorithm. (C-D) Network of top enriched human KEGG pathways for genes enriched from MSEA analysis.

KEGG orthology metabolic pathway analysis

Functional analysis was conducted with HUMAnN2²⁴ (version 2.8.1) by aligning clean reads (reads after removing host contamination) to KEGG protein database (version 2020.7). Furthermore, mapped KEGG genes were enriched to KO entries by MinPath²⁵. LEfSe analysis was performed to identify potentially enriched KEGG Orthology. We found that K07484 (transposase) and K07485 (transposase) were enriched in the MAFLD group (Fig. 5A,B). K06400 (site-specific DNA recombinase), and K19159 (antitoxin YefM) were enriched in the MAFLD with normal liver enzyme, but K07491 (REP-associated tyrosine transposase) was enriched in the MAFLD with abnormal liver enzyme (Fig. 5C,D).

KEGG orthology metabolic pathway analysis. (A-B) LDA scores and Heatmap showed significant pathway differences between the healthy and MAFLD groups. (C-D) LDA scores showed significant pathway differences between the MAFLD subgroups with or without liver enzyme abnormal.

Associations between metagenomics data and clinical parameters

To further explore the relationships between disturbances of gut microbiome and clinical variables, spearman correlation analysis was performed. We found that the differential bacterial microbiomes were generally associated with clinical variables (Fig. 6A). Furthermore, the representative microbiome and clinical factors with significant positive or negative correlations (P < 0.05, |Rho|> 0.4) between the indicated groups were shown by the correlation scatter plots as well (Fig. 6B). For example, Alistipes putredinis was significantly negatively correlated with GLU and GGT levels, Faecalibacterium prausnitzii was negatively correlated with TG, but Ruminococus gnavus was positively related to TG between the healthy controls, MAFLD subgroups with abnormal liver enzyme and normal liver enzyme, respectively. Collectively, these results suggested that Ruminococus gnavus might promote MAFLD progression. In contrast, Alistipes putredinis and Faecalibacterium prausnitzii may be useful in the symptomatic relief of MAFLD.

Relationship of gut microbiome and clinical parameters. (A) The association between species level bacteria (metagenomics) and clinical variable using spearman correlation. (B) The correlation of representative microbiome and clinical parameters in different groups. (C) The association between genus level bacteria (16S) and clinical variable using spearman correlation. (D) The correlation of representative microbiome and clinical parameters in different groups. *P < 0.05.

Associations between KEGG Orthology, genus by 16S sequencing and clinical parameters

Correlation analysis between the KEGG Orthology and clinical variables in both MAFLD and healthy group showed that TG was negatively associated with K00558, K02003 and so on (Fig. S1A,B).

On the other hand, gut microbiome genus with clinical factors had much different patterns of associations (P < 0.05, |Rho|> 0.4). As shown in Fig. 6C,D, Alistipes was negatively related to GLU, ALT, and GGT, UGG_002 was negatively related to ALP, TG, and GGT, Faecalibacterium was negatively related to TG.

Associations between gut microbiome species with the degree of abnormal liver enzymes

We re-analyzed the liver enzymes levels of these patients, GGT, ALT, and AST levels above the upper limit of normal (ULN) were included as three indexes for the abnormal liver enzymes subgroups. Four types of abnormal liver enzymes set as: 1N (1ULN, at least 1 index reached 1xULN AND under 1.5xULN), 1.5N (at least 1 index reached 1.5xULN AND under 2xULN), 2N (at least 1 index reached 2xULN AND under 3xULN), and 3N (at least 1 index reached 3xULN AND above).

The gut microbial species α‐diversity analysis revealed that the Shannon, Simpson and InvSimpson diversity indexes were lower in 1–1.5N group than in 2-3N, healthy, and normal liver enzymes groups; however, there was no significant difference between other groups (Fig. 7A). The microbiome compositions between the 1-1.5N and normal groups (PERMANOVA R² = 0.0142, P = 0.021), 2-3N and healthy group (PERMANOVA R² = 0.0335, P = 0.002) were significantly different (Fig. 7B).

The comparison of gut microbiome with the degree of abnormal liver enzymes. (A) Alpha diversity indices of the gut microbiome in different group in species level using metagenomics data. (B) The decomposition visualization of the beta diversity among the 4 groups in species level using metagenomics data. (C) Relative abundance of differential abundance microbiome species in different groups. (D) Alpha Diversity indices of the gut microbiome in different group in genus level using 16S data. (E) The decomposition visualization of beta diversity among the 4 groups. (F) Relative abundance of differential abundance microbiome genus in different groups. *P < 0.05, **P < 0.01, ***P < 0.001, ns P > 0.05.

To further investigate key microbiome related to the severity of MAFLD, we performed the relation analysis among microbiome genus and species with disease severity. According to the three subgroups above, including 1-1.5N, 2-3N, and normal liver enzyme group, we calculated the intersection of the three sets and obtain 3 species the 1-1.5N and 2-3N (Fig. 7C). Further analysis revealed that Odoribacter splanchnicus and Bacteroidales bacterium ph8 were significantly enriched within healthy group in contrast to the subgroups of 1-1.5N and 2-3N. Moreover, Fusobacterium mortiferum was enriched within the subgroups of 1-1.5N and 2-3N in contrast to healthy group.

Associations between gut microbiome genus with the degree of abnormal liver enzymes

The gut microbial genus α‐diversity analysis revealed that the Shannon, Simpson and InvSimpson diversity indexes were lower in 1-1.5N subgroup than in healthy groups; however, there was no significant difference between other groups (Fig. 7D). The microbiome composition difference between the MAFLD normal enzyme subgroup and healthy groups (PERMANOVA R² = 0.0188, P = 0.034), 1-1.5N and healthy groups (PERMANOVA R² = 0.0252, P = 0.004), 2-3N and healthy group (PERMANOVA R² = 0.0271 P = 0.039) were significantly different, indicating that there were composition differences between healthy and other groups (Fig. 7E). For gut microbiome genus, we obtained 5 genus from the 1-1.5N and 2-3N (Fig. 7F) subgroups. We found that Adlercreutzia, Alistipes, and Odoribacter were significantly enriched within healthy group in contrast to the 1-1.5N and 2-3N subgroups. Lachnospiraceae NK4A136 was enriched within healthy group in contrast to the 2-3N subgroup. Of note, Dorea was enriched in the 1-1.5N and 2-3N subgroups in contrast to healthy group.