We screened 1770 participants from the physical examination center of the First Affiliated Hospital of Zhengzhou University between March 2018 and May 2019. Inclusion criteria were as follows: (1) male over 18 years old; (2) clinically proven MAFLD according to the Asian Pacific Association for the Study of the Liver guidelines; and (3) legally competent to provide written informed consent. Exclusion criteria were: (1) the diagnosis of other liver diseases except MAFLD (for example hepatitis B or C, autoimmune hepatitis, hemochromatosis, or drug-induced hepatitis); (2) a history of inflammatory bowel disease; (3) treatment with antibiotics, probiotics or proton pump inhibitors within 2 months prior to the study; (4) the presence of an implanted orthopedic steel device or electronic medical device; or (5) the existence of a tumor. All patients provided written informed consent. The study was approved by the ethics committee from the First Affiliated Hospital of Zhengzhou University (2018-KY-56 and 2018-KY-90).

Our study was a retrospective cross-sectional study that was designed to evaluate the gut microbial species in MALFD patients with differing LSM values. A flow diagram describing participant selection is shown in Figure 1. After selection, there were 17 patients in the case group and 68 in the control group. The groups were matched by age, gender, waist circumference (WC), body mass index (BMI), dietary habits, and alcohol consumption. Patient characteristics are summarized in Table 1.

The subjects’ demographic information, behavioral risk factors, weight, height, WC, and blood pressure were collected by trained staff in the First Affiliated Hospital of Zhengzhou University. Samples to measure white blood cells (WBC), alanine transaminase (ALT), aspartate aminotransferase (AST), gamma-glutamyl transpeptidase (GGT), albumin (ALB), total bilirubin (TBIL), total cholesterol (TC), triglycerides (TG), high-density lipoprotein cholesterol (HDL), low-density lipoprotein cholesterol (LDL), fasting plasma glucose (FPG), glycated hemoglobin A1c (HbA1c), creatinine (Cre), and serum uric acid (SUA) were collected from the patients after 12 hours of fasting. On the same day, stool samples were obtained, immediately stored at -20°C, and frozen at -80°C within 30 minutes of collection. The samples were then tested.

A Toshiba Color Doppler Ultrasound System (APLI0500 TUS-A500) was used to scan patient livers. Abdominal ultrasound diagnosis of fatty liver required two of the following three observations: (1) the near-field echo of the liver was diffusely enhanced, and the echo was stronger than that of the kidney; (2) the structure of the intrahepatic duct was unclear; and (3) the far-field echo of the liver was gradually attenuated. Liver stiffness (LSM) and fat attenuation (CAP) were detected by FibroTouch (FT-C, Wuxi, China). The scans were performed by two experienced doctors.

Genomic DNA extraction from a total of 85 fecal samples was performed using the MagPure Stool DNA KF kit according to the manufacturer’s instructions. DNA libraries were constructed based on DNA nanospheres and shotgun metagenomic sequencing through binding probe anchor synthesis for all samples (MGI2000, MGI, Shenzhen, China). The overall accuracy (≥ 0.8) was used to evaluate the sequencing reads in order to filter out low-quality reads.

The intestinal microbiota was analyzed through classification, annotation, and quantification, based on MetaPhlAn2 with default settings (Truong et al., 2015). This analysis included bacteria, archaea, eukaryotic viruses, and fungi. The number of fungi present was small and is shown in Supplementary Figure 1. We then used HUMAnN2 (HMP Unified Metabolic Analysis Network 2) to further generate taxon-specific community function profiles (Franzosa et al., 2018). Neither archaea nor eukaryotic viruses were present in the final dataset; therefore, we did not include these in further analysis.

Statistical analyses were performed using the program R (version 4.0.2; www.r-project.org). Categorical variables are represented by counts, and chi-square tests were used for differential analyses. Continuous variables are expressed as means ± standard deviations or as medians (interquartile ranges). The Kolmogorov–Smirnov test was to test for normality, and Levene’s test was used to test the homogeneity of variances. Their results are presented in Supplementary Table 1. Depending on the results, either the Student’s t-test or the Mann–Whitney test were used. P < 0.05 was considered statistically significant. Permutational multivariate analysis of variance (PERMANOVA) and redundancy analysis (RDA) were performed to confirm whether liver stiffness was the most important influencing factor. The Shannon index and observed species number (obs) index of each sample was calculated using the vegan package in R. Principal coordinate analysis (PCoA) was performed using ade4 in R for visual analysis. Differences in the microbiota at the genus levels and pathways were analyzed using STAMP (version 2.1.3). Differences between the groups were calculated using the Benjamini–Hochberg false discovery rate (FDR). Species with low occurrences (occurrence rate < 10%) were removed before analyzing the microbiota. Spearman’s correlation analysis was used to analyze the correlations between the patients’ microbiota and other covariates, and the corrplot package was used for visualization.

A total of 1770 participants qualified for the study. Of these, 1607 either did not meet the inclusion criteria or failed to complete the liver scan with FibroTouch, and so were excluded. The remaining 163 participants completed the study and their LSM was assessed. Those with an abnormal LSM (LSM ≥ 7.4) were selected for the case group. To form a matched control group, 68 men were selected (all with LSM < 7.4 kPa) from the participants. The final cohort therefore included data from 85 patients (Table 1). A flow diagram depicting patient selection is shown in Figure 1A. The incidence of type 2 diabetes (T2DM) was higher in the case group than in the control group (P < 0.05). The median levels of alanine aminotransferase (ALT) and aspartate aminotransferase (AST) in the case group were 40.00 (interquartile range: 23.50–55.50) U/L and 26.00 (interquartile range: 23.00–33.50) U/L respectively. In both cases, the levels were higher than in the matched control group (P < 0.05). The liver stiffness (mean LSM = 8.94 ± 1.12 kPa) and fat attenuation index (median CAP = 300.30 dB/m; interquartile range: 265.98–326.05) were also significantly higher in the case group compared with the control group (P < 0.05). The other variables showed no statistically significant differences between the two groups.

PERMANOVA was used to analyze the basic attributes of the 85 participants (age, WC, systolic blood pressure (SBP), diastolic blood pressure (DBP), BMI, regular diet, dietary habits, consumption of whole grains, consumption of yogurt, smoking, alcohol consumption, exercise level, T2DM, and hypertension (HP)). In both univariate and multivariate analyses, a high LSM had the greatest impact on the gut microbiota of the study participants (P < 0.05; Table 2).

We constructed an RDA diagram to examine the relationship between the microbiota and participants’ dietary habits and individual attributes (Figure 1B).

We next analyzed the bacteria present in the gut microbiota. The case group was dominated by the phyla Bacteroidetes, Firmicutes and Actinobacteria. In the control group, Bacteroidetes and Firmicutes were dominant. We observed seven bacterial species that differed between the case group and the control group. Furthermore, we observed two fungi (Naumovozyma unclassified and Saccharomyces cerevisiae) present only in the control group (Figure 1C). The bacteria Prevotella copri, Phascolarctobacterium succinatutens, Eubacterium biforme, and Collinsella aerofaciens were all increased in the case group compared with in the control group (P < 0.05). Conversely, Bacteroides coprocola, Bacteroides stercoris and Clostridiales bacterium 1_7_47FAA were all reduced in the case group (P < 0.05). We then assessed the α-diversity of gut microbiota using the Shannon diversity index and the observed species number (obs) index, both of which showed no statistically significant difference between the two groups (P > 0.05). To look at β-diversity, we performed principal coordinate analysis (PCoA1) using the Bray–Curtis and the Pearson distance tests. These analyses determined that there was no statistically significant difference between the case and control groups (P > 0.05). However, the composition of the gut microbiota from patients in the case group clustered independently compared with patients in the control group (Figure 2).

Next, we used Spearman’s correlation analysis to examine the correlation between serological indicators and the abundances of certain species within the gut microbiota. The level of Phascolarctobacterium succinatutens was positively correlated with LSM (P < 0.05). Eubacterium biforme was also positively correlated with LSM (P < 0.05), but was inversely correlated with levels of high-density lipoprotein cholesterol (HDL; P < 0.01). Clostridiales bacterium 1_7_47FAA was negatively correlated with white blood cell count (WBC; P < 0.05) but was positively correlated with fasting plasma glucose levels (FPG; P < 0.05). Collinsella aerofaciens was positively correlated with both LSM and WBC (P < 0.05) yet was negatively correlated with HDL (P < 0.05). Bacteroides stercoris was inversely correlated with low-density lipoprotein cholesterol (LDL) and LSM (P < 0.05) but positively correlated with WC, glycated hemoglobin A1c (HbA1c; P < 0.05) and CAP (P < 0.01). The fungus Naumovozyma unclassified was positively correlated with SBP, albumin levels (ALB; P < 0.01), and WBC (P < 0.05). The fungus Saccharomyces cerevisiae was positively correlated with WC and BMI (P < 0.05) (Figure 3).

We used the MetaCyc database to construct functional profiles for each sample. After removing low abundance pathways, a total of 32 pathways were found to be significantly different between the two groups (Figure 4; Supplementary Table 2). Of these pathways, 16 were enriched in the case group compared with the control group. Among these, 14 pathways (ARO-PWY, PWY-6163, PWY-5097, PWY-6151, PWY-2942, POLYAMINSYN3-PWY, PWY-6897, PWY-7357, PWY-6700, PWY-5686, PWY-7221, COMPLETE-ARO-PWY, PWY-7039, and PWY-7219) were related to biosynthesis of key metabolites, including amino acids, vitamins, nucleosides, and nucleotides. The remaining two, PWY-7117 and PWY-241, were related to photosynthesis (Figure 4).

We then applied Spearman’s correlation analysis to analyze the correlation between bacterial species and microbial pathways (Figure 5). The level of Clostridiales bacterium 1_7_47FAA was positively correlated with PWY-5154, PYRIDOXSYN−PWY, PWY0−845, PWY−5101 (P < 0.05), and ARGININE−SYN4−PWY (P < 0.01). Bacteroides stercoris was positively correlated with PWY0−845, ARGININE−SYN4−PWY, PWY−5101 (P < 0.001), and PYRIDOXSYN−PWY (P < 0.01), but negatively correlated with PWY−7357 (P < 0.05). Collinsella aerofaciens was positively correlated with PWY−241, PWY−7117 (P < 0.001), and POLYAMINSYN3−PWY (P < 0.01), yet negatively correlated with PYRIDOXSYN−PWY, PWY0−845 and, ARGININE−SYN4−PWY (P < 0.05). Phascolarctobacterium succinatutens was positively correlated with POLYAMINSYN3−PWY (P < 0.05). Bacteroides coprocola was negatively correlated with REDCITCYC, PWY−5154 (P < 0.001), HEXITOLDEGSUPER−PWY, ANAEROFRUCAT−PWY, and PWY−7039 (P < 0.05). Eubacterium biforme was positively correlated with POLYAMINSYN3−PWY (P < 0.001), but negatively correlated with NAGLIPASYN−PWY, ANAEROFRUCAT−PWY, PWY−384, and PYRIDOXSYN−PWY (P < 0.05). Prevotella copri was negatively correlated with PWY−5695 (P < 0.001), COLANSYN−PWY, HEXITOLDEGSUPER−PWY, PWY−3841, 1CMET2−PWY, PWY−6703, and PWY−5101 (P < 0.05). The fungus Naumovozyma unclassified was positively correlated with PWY−7039 (P < 0.05). The fungus Saccharomyces cerevisiae was negatively correlated with PWY−241 (P < 0.05).

The correlation between MetaCyc pathways that showed significant differences between the case and control groups and other clinical features was then examined using Spearman’s correlation analysis. The MetaCyc pathways were clustered based on MetaCyc pathway hierarchy and constructed into a heatmap (Figure 6). We observed that some pathways were negatively correlated with LSM in the control group, while the majority of pathways in the case group, as well as a subset of pathways in the control group, were negatively correlated with both HbA1c and age.

Finally, we used a random forest model to build a predictive model for liver fibrosis (Figure 7). The model was constructed using seven bacterial species (Prevotella copri, Phascolarctobacterium succinatutens, Eubacterium biforme, Collinsella aerofaciens, Bacteroides coprocola, Bacteroides stercoris and Clostridiales bacterium 1_7_47FAA). The area under the curve (AUC) for the model was 0.728 (95%; confidence interval: 0.574–0.881). The model showed good diagnostic performance.