Role of intestinal flora in the development of nonalcoholic fatty liver disease in children

A total of 114 participants were recruited for the study between January 2019 and December 2021, comprising 79 patients with NAFLD and 35 healthy controls (HCs). The participants were aged 7–18 years, with information on various parameters such as sex, birth date, height, and weight collected from their professionally trained parents. Recruited in the NAFLD group were patients whose body mass index (BMI) was higher than the 95th percentile and defined through a clinical evaluation and liver ultrasonography (10). The BMI values were used to classify overweight and obesity among participants, with BMI calculated as the ratio of weight (kg) to height squared (m²). Participants with a recent history of antibiotic intake within 2 months or probiotics use within a month before enrollment were excluded. The liver function tests and lipid profile assessment were conducted on fasting blood samples, and the fecal samples were also collected (11). The diagnostic criteria for metabolic disorders of NAFLD were as follows (12): (i) central obesity or overweight; (ii) antidiabetic treatment or fasting glucose level ≥ 7.0 mmol/L; (iii) antihypertensive treatment or elevated collected blood pressure with ≥95th percentile in children of the same age and sex; (iv) treatment with lipid-lowering agents or triacylglycerols ≥ 1.47 mmol/L; and (v) treatment with lipid-lowering agents or high-density lipoprotein cholesterol < 1.1 mmol/L in men or <1.3 mmol/L in women. Fifteen children with NAFLD were reassessed for their BMI, and fecal samples for metagenomic sequencing were obtained following 3–6 months of nonmedical treatment such as diet control and increased exercise.

The liver biopsies were collected from children whose parents or guardians provided consent for liver puncture to assess NAFLD activity score (NAS) and liver fibrosis stage. The liver tissues of 13 children were obtained through ultrasound-guided liver biopsy and fixed in formaldehyde. The specimens were embedded in paraffin, after which serial sectioning was performed, and specimens were stained with hematoxylin-eosin (HE) and reticular fiber and/or Masson trichrome staining. The slides were then examined under a light microscope. Two pathologists scored the pathological diagnosis results at the same time; in case of disagreement, a senior physician was invited. NAFL was assessed according to NAS diagnosis criteria of NAFLD on pathological liver tissue (or liver biopsy) (13). NAFL was defined as liver cell steatosis without any indication of steatohepatitis, with or without liver fibrosis; NASH was defined as the presence of hepatic steatosis with necroinflammation and hepatocellular injury with or without fibrosis. NAFLD with fibrosis can refer to either NAFL or NASH associated with fibrosis in port, periportal, perisinusoidal, or portal-portal bridging fibrosis (13). The diagnostic criteria for NAFLD and NASH, including the scoring rules for NAS, the stage of liver fibrosis, and the pathological diagnosis of NASH, refer to the Expert Consensus on the Diagnosis and Treatment of NAFLD in Children and the study by Kleiner et al. (13, 14).

Out of the 79 children diagnosed with NAFLD, our analysis included only 13 children for whom parental consent was obtained for liver puncture. Following the evaluation of the NAS, we identified eight cases of NASH and five cases of NAFL. As shown in Fig. 1, the NAS was created as an unweighted score for steatosis (0–3), lobular inflammation (0–3), and ballooning (0–2). A four-scale grading (from 0 to 3) was used for grading steatosis: normal liver (grade 0) contains fat in <5% of hepatocytes, while grade 1 steatosis refers to <33% steatotic hepatocytes. In grade 2 and 3 steatoses, fat is present in at least 33% or 66% of hepatocytes, respectively. The grading system of lobular inflammation was determined by the number of inflammatory foci in a given area, from grade 0 to 3 (0–3 points, none, <2, 2–4, and >4, respectively). The grading of ballooned cells was based mainly on number (0–2 points, none, few, and many, respectively). The total NAS score ranged from 0 to 8, where scores of 0–2 can exclude NASH, 3–4 suggested possible NASH, and scores of 5–8 were used to diagnose NASH. Additionally, the progression of fibrosis was described as follows (Fig. 1): zone 3 perivenular and perisinusoidal represented stage 1 (with a subdivision into 1a and 1b according to the amount of the deposit), and stage 1c was isolated periportal fibrosis; stage 2 included portal and central fibrosis without bridging fibrosis; stage 3 was bridging fibrosis; and stage 4 was cirrhosis.

After overnight fasting, 4 mL of peripheral blood was drawn from the cubital vein of each participant (15). The blood samples were centrifuged at 2,500 g for 10 min using an analyzer (KHB, Shanghai, China). The resulting supernatants were collected and stored for analysis. The serum levels of fasting blood glucose, uric acid, thyroid-stimulating hormone (TSH), alanine aminotransferase (ALT), and aspartate aminotransferase (AST) were measured using a Beckman Coulter AU5800 analyzer (CA, USA).

The fecal specimens from participants were collected in sterilized 2 mL tubes containing a protective solution and frozen at −80°C until DNA extraction. Total DNA was extracted from the fecal samples (250 mg, wet weight) using a PSP Spin StooL DNA Kit/PSP Spin Stool DNA Plus Kit (Stratec Molecular GmbH, Berlin, Germany) following the manufacturer’s protocols (10).

The stool samples were processed at MyGenostics Co. Ltd. (Beijing, China). The 16S rRNA amplicon sequencing of the V3–V4 gene region was used on an Illumina MiSeq platform to detect and quantify bacterial taxa in the fecal samples of HCs and children with NAFLD. The paired-end 16S sequences were merged and filtered for quality using the laser flash method. The Quantitative Insights into Microbial Ecology (QIIME, version 2.0) software suite was used to analyze the sequences, and the QIIME tutorial (http://qiime.org/↗) was followed with some modifications. The sequences were compared with the 97% reference data set of the SILVA database (SILVA project, version 138.1, http://www.arb-silva.de/↗). Those that did not match any reference sequences were clustered into de novo operational taxonomic units (OTUs) at 97% similarity using the UCLUST algorithm. This allowed us to accurately assign taxonomical classifications to our sequencing reads and perform downstream analyses with confidence. To obtain the population profiles of common and abundant taxa, we constructed a distance matrix based on the abundance of OTUs. Hierarchical clustering was performed using the unweighted pair group method with arithmetic mean (average linkage clustering) on this distance matrix. A Newick tree format was obtained using the QIIME package. Furthermore, QIIME was used to analyze the alpha diversity (observed species, Chao1, Shannon, and Simpson) and beta diversity [principal coordinate analysis (PCoA)]. Linear discriminant analysis effect size (LEfSe) was applied to compare significant differences in taxa between patients with NAFLD and HCs using a P value < 0.05.

In this study, we utilized the random forest supervised learning algorithm on 16S data sets to distinguish between HCs and NAFLD samples. We set the ntree parameter to 500 for the number of trees to construct. Two-thirds of all data were utilized as the training set, while the remaining third constituted the test set. To reduce the possibility of confounding variables due to correlated data, we normalized the filtered relative abundances. The model was trained using 114 stool samples, 35 of which were HC and 79 were NAFLD samples. The class error calculated using the model was 0.3714. The model accurately identified 22 HC samples. For the NAFLD samples, the model could predict all samples precisely. The estimated error rate of the model was 0.11404, while the baseline error for random guessing sat at 0.30702. The ratio of baseline error to observed error for our modeling was 2.69231. We used the R programming language to calculate both area under curve (AUC) and the top 10 key indicators generated by the model.

To investigate the bacterial species spectrum associated with NAFLD clinical diagnosis, we conducted metagenomic sequencing on 79 fecal samples obtained from children with NAFLD. The stool specimens were submitted to MyGenostics Co. Ltd. (Beijing, China) for analysis. During metagenomic sequencing, 20 µg of microbial metagenomic DNA was extracted and purified. Also, a sequencing library was prepared using a NEXTflex Rapid DNA Sequencing Kit (NOVA-5144-02, Bioo Scientific, TX, USA) following the manufacturer’s protocols. After library preparation, Qubit 2.0 was used to preliminarily quantify the library, which was then diluted to 2 ng/µL. Agilent 2100 was used to detect the insert size of the library, and the quantitative polymerase chain reaction method was used to accurately measure the effective concentration of the library (>3 mmol/L) to ensure the quality of the library. The constructed library was sequenced on an Illumina NovaSeq 6000 platform. MOCAT preprocessing was employed to ensure the quality of data, which involved filtering out sequences with low quality, insufficient length, or containing errors. The resulting clean data from each sample were then used for metagenomic assembly. Subsequently, all the samples were assembled, with the objective of generating reads longer than 300 bp for further analysis. The genes were combined and clustered using CD-Hit software, ultimately resulting in 2,053,172 nonredundant genes.

The gene catalog was used to compare the sequence with the MicroNR library to obtain species annotation information of unigenes employing bowtie. Subsequently, the gene abundance table was merged to generate a species abundance table with different classification levels. The overall distinctions in microbial composition were evaluated by PCoA analysis of the Bray-Curtis distance. LEfSe was performed, and a significance level of P < 0.05 was used to detect significant variations in taxa between patients with NAFL and NASH. The sequences obtained from metagenomic sequencing and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database were compared using BLAST (version 2.2.28+) to identify the metabolic functions of gut microbiota within the samples. The predicted KEGG Orthology-Based Annotation System 2.0 was used for this purpose. The relationship between the structure of intestinal microbiota and NAFLD-related indexes was analyzed using Spearman’s correlation.

The present study used SPSS 22.0 software and an R environment during statistical analysis. Before conducting the statistical analyses, the normality was evaluated for all numerical variables per previous recommendations (16). Generalized linear models were employed to compare continuous data across groups. The chi-square analyses were conducted to compare frequencies among groups, and the independent-sample t-tests were used to compare age and BMI across groups. The Kruskal-Wallis test was used to investigate the differences in taxonomic composition and function annotations between groups. Alpha diversity was calculated using the Shannon index. In contrast, differential diversity analyses were conducted using either the Wilcoxon rank-sum test or the Kruskal-Wallis test, depending on the number of groups. Kruskal-Wallis H-tests and Welch’s t-tests were used to calculate P values, and false discovery rates were employed to control for multiple hypothesis testing. Mann-Whitney U tests were used to compare the abundance of each phylum and genus among groups. For phyla and genera in which Student’s t-tests were employed, the abundance values were ranked and used in regression models. Additionally, each multiple linear regression model was checked for normality using predicted probability plots, and the absence of multicollinearity was assessed using variance inflation factor values.

This study comprised 114 children between 7 and 18 years old, among which 79 were with NAFLD and 35 were HCs. Height, weight, and BMI were measured, and no significant differences were observed in sex or age between the two groups. However, the BMI of patients with NAFLD was significantly higher than that of HCs. The anthropometric and demographic data are presented in Table 2.

The gut microbiomes of 79 patients with NAFLD and 35 HCs were analyzed using 16S rRNA sequencing to explore the differences in intestinal flora diversity and the abundance of various taxa. A total of 16,570,490 sequencing reads were obtained from the 114 samples. The fecal microbiotas observed in the samples were primarily composed of the phyla Firmicutes, Bacteroidetes, and Proteobacteria (Fig. 2A). Only the abundance of Proteobacteria was found to increase in pediatric patients with NAFLD compared with HCs. The dominant genera observed in both groups were Bacteroides, Faecalibacterium, Prevotella 9, and Bifidobacterium (Fig. 2B). However, the abundance of Bacteroides and Prevotella 9 both increased in pediatric patients with NAFLD, whereas the abundance of Faecalibacterium decreased. Differences in the abundance of several bacteria were observed between NAFLD and HC groups, and these bacteria were planned for further analysis. No significant differences in species richness and evenness were observed between HCs and patients with NAFLD in median alpha diversity (observed species, Chao1, Shannon, and Simpson diversity indices) (Fig. 2C). However, the observed species in patients with NAFLD were more dispersed than those observed in HCs. A PCoA analysis of the gut microbiota samples of the children is presented in Fig. 2D. The statistical models revealed that the composition and abundance of flora species in patients with NAFLD and HCs were similar.

The analysis of bacterial abundance among groups revealed a significant difference in the abundance of several species of bacteria (Fig. 3). Specifically, the abundance of 24 species of bacteria, including Lachnospiraceae, Proteobacteria, Ruminococcus gnavus group, Fusobacteria, Fusobacteriales, Lachnoclostridium, Enterobacteriaceae, Enterobacteriales, and Escherichia-Shigella, increased in the gut microbiota of pediatric patients with NAFLD. Conversely, the abundance of 51 species of bacteria, including Ruminococcaceae, Faecalibacterium, Bifidobacteriaceae, Bifidobacteriales, Prevotella 9, Prevotellaceae, Ruminococcus 1, and Veillonellaceae, was significantly lower in patients with NAFLD compared with HCs (Fig. 3). Additionally, Methanobacteria belonging to Euryarchaeota was a characteristic genus in HCs. Conversely, Fusobacteria was the characteristic genus in patients with NAFLD.

Eight patients with NASH and five with NAFL were identified among children who consented to liver puncture, and the pathological findings were consistent with the diagnosis of NAFLD. Patient BMI, fasting blood glucose, uric acid, TSH, and two liver enzyme indices were measured and compared between the NAFL and NASH groups. Table 3 illustrates the anthropometric and clinical features of the patients. Patients with NASH had higher levels of AST and TSH than those with NAFL.

The samples from patients with NAFL and NASH underwent metagenomic sequencing to determine the diversity and significant differences in intestinal flora between the two groups. We discovered that the number and evenness of flora species, as well as the Shannon and Simpson diversity indices, decreased in patients with NAFL compared with NASH (Fig. 4A). The non-metric multidimensional scaling (NMDS) plot of the gut microbiota samples revealed that the abundance and composition of flora species in patients with NASH were similar to those in patients with NAFL (Fig. 4B). Further analysis using LEfSe identified several significantly different genera and species between the two groups. Patients with NASH had a higher abundance of Alistipes, whereas Peptostreptococcaceae noname had a lower abundance at the genus level (P < 0.05, Fig. 4C).↗ The abundance of Bacteroides uniformis, Lachnospiraceae bacterium 7_1_58FAA, Eubacterium ventriosum, and unclassified Roseburia increased significantly in patients with NASH compared with those with NAFL at the species level (P < 0.05, Fig. 4D).

We conducted metagenomic sequencing on all 79 NAFLD samples to analyze the correlation between NAFLD-related indices and microbiota in children. The abundance of certain species positively correlated with NAFLD-related indices. The abundance of Streptococcus mitis, Streptococcus oralis, Streptococcus pneumoniae, Gemella unclassified, Collinsella unclassified, Granulicatella adiacens, and Streptococcus cristatus positively correlated with BMI (Fig. 5). GLU positively correlated with the abundance of Weissella confusa, Lachnospiraceae bacterium 2_1_58FAA, Roseburia intestinalis, Lachnospiraceae bacterium 6_1_63FAA, and Paraprevotella xylaniphila (Fig. 5). ALT positively correlated with the abundance of Eubacterium ruminantium group, Ruminococcaceae NK4A214 group, Erysipelotrichaceae UCB-003, Ruminococcus inulinivorans, and Coprococcus comes. AST positively correlated with the abundance of Barnesiella intestinihominis, Eubacterium redale, R. inulinivorans, and Akkermansia muciniphila (Fig. 5). Only the abundance of unclassified Neisseria positively correlated with TSH (Fig. 5). Uric acid (UA) negatively correlated with the abundance of Ruminiclostridium 9, Collinsella aerofaciens, and Eubacterium eligens, whereas Collinsella intestinalis, Collinsella unclassified, Bacteroidetes xylanisolvens, and Megamonas rupellensis showed a positive correlation (Fig. 5). These findings suggested a correlation between gut microbiota dysbiosis and NAFLD development, with these bacteria being considered NAFLD-related taxa.

We conducted a further examination of disruptions in gut microflora associated with specific metabolic pathways in NAFLD. These pathways included glucose metabolism, fatty acid metabolism, short-chain fatty acid (SCFA) metabolism, and folate pathways. Utilizing metagenomic sequencing data from 79 NAFLD patients, we assessed the correlation between bacterial taxa and these metabolic pathways. As shown in Fig. 6, the abundance of Turicibacter sanguinis, Eubacterium dolichum, and Bacteroides clarus was positively correlated with xylose degradation IV (PWY-7294). The abundance of Weissella confusa, Actinomyces odontolyticus, S. mitis, S. oralis, S. pneumoniae, and Streptococcus australis was positively correlated with N10-formyl-tetrahydrofolate biosynthesis and folate transformations II (PWY-3841), whereas negatively correlated with the abundance of Clostridiales noname. The abundance of Citrobacter freundii and Citrobacter unclassified was positively correlated with fatty acid beta-oxidation V (PWY-6837), fatty acid salvage (PWY-7094), and L-glutamate degradation IV (PWY-4321). The abundance of Neisseria subflava, Neisseria flavescens, Fusobacterium periodonticum, Lautropia mirabilis, and T. sanguinis showed a positive correlation with the superpathway of taurine degradation (PWY-1541).

The abundance of different bacterial taxa was strongly associated with the metabolism of SCFAs (Fig. 6). Bacteria such as Clostridiales bacterium 1_7_47FAA, Escherichia coli, unclassified Escherichia, Veillonella parvula, R. gnavus, and C. freundii were found to be strongly linked to the tricarboxylic acid (TCA) cycle VII, which was involved in the production of acetate. Bacteria such as unclassified Oscillibacter, Dorea longicatena, Adlercreutzia equolifaciens, Holdemania filiformis, Atopobium, Clostridium leptum, Roseburia hominis, R. obeum, E. dolichum, Odoribacter splanchnicus, Bacteroides intestinalis, and Veillonella dispar were strongly linked to pyruvate fermentation to acetate and lactate II (PWY-5100). The abundance of Enterobacter cloacae and B. clarus was found to be strongly associated with acetyl-CoA fermentation to butyrate. Similarly, the abundance of Solobacterium moorei, unclassified Granulicatella, and Clostridiales bacterium 1_7_47FAA was strongly linked to the hexitol fermentation to lactate, formate, ethanol, and acetate (P461-PWY). These findings indicate a close relationship between the composition of fecal bacteria and specific metabolic pathways within the gut microbiota. Perturbations in bacterial composition could potentially impact the development of NAFLD through microbiota metabolism.

We developed a stochastic forest machine algorithm to construct a prediction model that could discriminate between patients with NAFLD and HCs based on the identified bacterial taxa. We performed 16S rRNA sequencing analysis of the gut microbiomes of patients with NAFLD and HCs. Ten bacterial taxa, including Sneathia, Lachnospiraceae ND3007 group, Prevotella, Lachnoclostridium, Escherichia-Shigella, Prevotellaceae UCG-001, Haemophilus, Lactobacillus, Ruminiclostridium 6, and Coprococcus 2, were selected to predict the risk of NAFLD in children (Fig. 7A). The prediction accuracy of the model was assessed using receiver operating characteristic (ROC) curves, which are presented in Fig. 7B. The model achieved an AUC value of 0.969, with a sensitivity of 0.914 and a specificity of 0.949.

In this study, we examined the alterations in the diversity of the gut microbiota and abundance of NAFLD-related bacteria before and after nonpharmacological treatment to establish the correlation between gut microbiota and BMI reduction. Fifteen children diagnosed with NAFLD agreed to participate in the study. Five children recorded a decrease in BMI greater than 5% after the nonpharmacological intervention, and their treatment was labeled effective. Conversely, two children had a BMI increase greater than 5%, indicating ineffective treatment. The remaining eight children showed a decrease in BMI of less than 5%, which was deemed irrelevant in this study. Subsequently, we conducted metagenomic sequencing and the resultant analysis of microbial taxa on stool samples of children whose BMI increased or decreased by more than 5%. Figure 8A displays the study results, revealing an increase in the abundance of gut flora species in children with reduced BMI after nonpharmacological treatment. In contrast, those with increased BMI experienced a decrease. The abundance of Bifidobacterium longum and E. ventriosum increased significantly with nonpharmacological treatment (Fig. 8B). Furthermore, the abundance of Alcaligenaceae, Erysipelotrichaceae, Faecalibacterium prausnitzii, Subdoligranulum unclassified, Lachnoclostridium, R. intestinalis, Lachnospiraceae UCG-008, and C. comes all demonstrated significant increases following nonpharmacological treatment (Fig. 8B). Conversely, the study noted a decrease in the abundance of Haemophilus, E. coli, Escherichia unclassified, Ruminococcus torques, and Lachnospiraceae bacterium 2_1_58FA in patients with NAFLD with reduced BMI after treatment (Fig. 8C). These results suggested a correlation between the diversity of gut microbiota, the abundance of key species, and the reduction of BMI in NAFLD after nonpharmacological treatment.

Sequences analyzed in this study are available in GenBank with the accession number PRJNA914785↗.