A randomized controlled trial for response of microbiome network to exercise and diet intervention in patients with nonalcoholic fatty liver disease

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ A randomized controlled trial for response of microbiome network to exercise and diet intervention in patients with nonalcoholic fatty liver disease © The Author(s) 2022 Published version Cheng, Runtan; Wang, Lu; Le, Shenglong; Yang, Yifan; Zhao, Can; Zhang, Xiangqi; Yang, Xin; Xu, Ting; Xu, Leiting; Wiklund, Petri; Ge, Jun; Lu, Dajiang; Zhang, Chenhong; Chen, Luonan; Cheng, Sulin Cheng, R., Wang, L., Le, S., Yang, Y., Zhao, C., Zhang, X., Yang, X., Xu, T., Xu, L., Wiklund, P., Ge, J., Lu, D., Zhang, C., Chen, L., & Cheng, S. (2022). A randomized controlled trial for response of microbiome network to exercise and diet intervention in patients with nonalcoholic fatty liver disease. Nature Communications, 13, Article 2555. https://doi.org/10.1038/s41467-022-299680 2022 ARTICLE A randomized controlled trial for response of microbiome network to exercise and diet intervention in patients with nonalcoholic fatty liver disease Runtan Cheng1,2,3,13, Lu Wang4,13, Shenglong Le1,5, Yifan Yang 1,3, Can Zhao6, Xiangqi Zhang1,3, Xin Yang 3, Ting Xu3, Leiting Xu1,7, Petri Wiklund1,2, Jun Ge1,8, Dajiang Lu2,9, Chenhong Zhang 3,10✉, Luonan Chen4,11,12✉& Sulin Cheng 1,2,5✉ Exercise and diet are treatments for nonalcoholic fatty liver disease (NAFLD) and prediabetes, however, how exercise and diet interventions impact gut microbiota in patients is incompletely understood. We previously reported a 8.6-month, four-arm (Aerobic exercise, n=29; Diet, n=28; Aerobic exercise +Diet, n=29; No intervention, n=29) randomized, singe blinded (for researchers), and controlled intervention in patients with NAFLD and prediabetes to assess the effect of interventions on the primary outcomes of liver fat content and glucose metabolism. Here we report the third primary outcome of the trial—gut microbiota composition—in participants who completed the trial (22 in Aerobic exercise, 22 in Diet, 23 in Aerobic exercise +Diet, 18 in No Intervention). We show that combined aerobic exercise and diet intervention are associated with diversified and stabilized keystone taxa, while exercise and diet interventions alone increase network connectivity and robustness between taxa. No adverse effects were observed with the interventions. In addition, in exploratory ad-hoc analyses we find that not all subjects responded to the intervention in a similar manner, when using differentially altered gut microbe amplicon sequence variants abundance to classify the responders and low/non-responders. A personalized gut microbial network at baseline could predict the individual responses in liver fat to exercise intervention. Our findings suggest an avenue for developing personalized intervention strategies for treatment of NAFLD based on host-gut microbiome ecosystem interactions, however, future studies with large sample size are needed to validate these discoveries. The Trial Registration Number is ISRCTN 42622771. https://doi.org/10.1038/s41467-022-29968-0 OPEN A full list of author affiliations appears at the end of the paper. NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications 1 1234567890():,; Nonalcoholic fatty liver disease (NAFLD) is the most common chronic liver disease with prevalence estimates ranging from 25% to 45% in worldwide1,2. NAFLD is closely associated with type 2 diabetes (T2D)3, most likely via a common pathophysiological mechanism - insulin resistance4. Therefore, interventions targeting the synergistically pathogenic mediator of NAFLD and diabetes are likely to be a rational approach in the prevention and treatment of the comorbidity condition. The gut microbiota, which contributes to the metabolic health of the human host, may work as one of the targets for the treatment of NAFLD. Emerging data demonstrate that dysbiosis of the gut microbiota is linked with NAFLD and T2D5,6. Although the exact mechanism(s) require further elucidation, inflammation, damage to the intestinal membrane, and translocation of bacteria have all been suggested7. Drugs with antifibrotic or anti-inflammatory treatments for various stages of NAFLD in human trials is pending8. Currently, increased physical activity and dietary modifications are the only effective therapeutic options for NAFLD management, and the mechanisms of these interventions have been associated with modulation of gut microbiota and its metabolites9. For example, a low-carbohydrate diet (LCD) has been reported to improve fatty liver metabolism and to promote rapid shifts of the gut microbiota composition in NAFLD patients10. Regular aerobic exercise has been shown to reduce hepatic fat in obesity11,12, and some studies have shown that exercise increases gut microbial diversity13 and alters the composition and functional capacity of gut microbiota14,15. However, the underlying mechanism of the beneficial effects of exercise and diet on NAFLD, as well as the related metabolic disorders through regulation of the gut microbiota, still remain to be elucidated. One of the challenges in constructing the relationship between the gut microbiota and improvement of NAFLD during exercise and/or dietary interventions is that not all patients respond to the intervention in a similar manner. Earlier studies have shown that there was an individual variation in response to exercise, with some subjects experiencing greater improvement than others16. The low/non-responders could be up to 50% in terms of change of gut microbiota after exercise intervention17. Some researchers have suggested that the variation in responsiveness might be due to genotypic and phenotypic factors18,19. Thus, large inter-individual variance may mask the changes of the microbiome in studies using traditional statistical analytical methods. Currently, it is widely accepted that, in the microbial community, the keystone taxa are drivers of microbiome structure and function, and in particular, their interaction network, which plays an important role in microbial functions and disease progression20. Therefore, to establish an effective intervention strategy, it may be worthwhile to analyze the microbiota correlation network at the microbial community level, while assessing the physiological mechanisms at the individual level and further exploring individual microbial networks that underlie the differences between responders and low/nonresponders of various interventions. In the present study, we analyzed the composition and metabolic pathways of the gut microbiota, constructed a co-occurrence network at the population level, and developed personalized gut microbial networks. This individual network enabled us to identify the microbial signature and interaction of taxa at a personal resolution, and to further differentiate responders from low/ non-responders after intervention. Taken together, our study with relatively long-term intervention provides both statistical and sample-specific network insights into the complex gut microbial ecosystem in patients with NAFLD and glucose metabolism impairment. Results Participant characteristics and gut microbiota research design. This study was an 8.6-month, four-arm, randomized trial (Fig. 1a, Fig. S1) and participant characteristics have been reported in our previous publication21.Briefly, 115 participants were recruited from 7 health clinical service centers in the Shanghai Yangpu district. They were randomized into four groups: aerobic exercise intervention (AEX, n=29), fiber-enriched low-carbohydrate diet intervention (Diet, n=28), aerobic exercise combined with diet intervention (AED, n=29) and no intervention without guided exercise and dietary intake (NI, n=29). Of all participants, 85 individuals completed the intervention trial. For the primary outcome of liver fat content, we previously found that hepatic fat content (HFC) was significantly reduced in the exercise AEx (–24.4%), diet (–23.2%), and AED (–47.9%) groups by contrast to the 20.9% increase in the NI group (p< 0.001 for all) after intervention21. Of note, 91% of the subjects in the AED group decreased their HFC, and the corresponding figures were 68% in the AEx group, and 86% in the Diet group. In contrast, 72% in the NI group increased their HFC during the intervention period. However, for the primary outcome of glucose metabolism, no significant remission or progression of prediabetes was found between the intervention and NI groups based on the glycated hemoglobin A1c (HbA1c)21. These results indicated that our intervention was effective mainly for HFC reduction. We found that not all the subjects responded to the intervention in a similar manner. Therefore, we stratified the participants into responders (HFC decreased more than 5%) and low/nonresponders (HFC decrease less than 5% or increased) (Fig. 1b). For the primary outcome of gut microbiota composition, seventy-six subjectsprovidedpairedstoolsamplesatbaselineandafterthe intervention for 16 s rRNA gene sequencing. The demographics and clinical variables of those who have had gut microbiota results are presented to Table S1. In addition, to better understand the intervention impact on the function of microbiome, we selected a subset of the cohort (n=42) with the best or worst response, in terms of HFC reduction, to the intervention, and analyzed gut microbiota by the metagenomics data (Fig. 1c). The microbiota data and clinical parameters were used to (1) characterize the change of gut microbiota composition in response to interventions; (2) determine associations of the microbiome with the level of physical fitness, fat mass, serum biomarkers and short-chain fatty acids (SCFAs); (3) analyze the metabolic function shift of the microbiome in response to interventions; (4) discover the co-occurrence network features of the microbiome ecosystem before and after intervention; and (5) predict personalized response to intervention based on the baseline gut microbiota composition and network. Changes of structure and function of gut microbiota induced by exercise and/or dietary intervention of participants with comorbidity of NAFLD and prediabetes.Wefirst performed 16 S rRNA gene sequencing of fecal samples collected before and after interventions and obtained a total of 5421 amplicon sequence variants (ASVs). We found that the alpha diversity (Shannon index) of gut microbiota was significantly decreased in the NI group (p=0.045) with large individual variance by contrast to the intervention groups, which maintained their diversity after the intervention (Fig. 2a, NI vs. AED p=0.011, NI vs. AEx p=0.007 and NI vs. Diet p=0.025, respectively, analysis of variance with repeated measures adjusted for change of body weight, baseline value and intervention duration). The changes in microbial diversity may reflect the natural progression of NAFLD in the NI group. Moreover, the principal coordinates analysis (PCoA) based on weighted UniFrac distances shows the ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29968-0 2NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications significant difference between NI and the other groups after the interventions (p< 0.05, PERMANOVA, Fig. 2b). These results indicate that microbial diversity deteriorates with the increased HFC, while exercise and dieting may help to maintain the diversity of the gut microbiota. To avoid only included the relative abundance data in the analysis which may introduce the bias, we then measured the absolute abundance of microbiome by quantitative real-time PCR (qPCR) targeting the 16 S rRNA gene. We found that there were no significant differences between baseline and after intervention as well as among the groups (Fig. S2a). Moreover, the total bacterial content of each subject in the same group at two time points also did not differ significantly (Fig. S2b–e). Furthermore, we performed Linear discriminant analysis Effect Size (LEfSe) analysis22 of ASVs that appeared in more than 20% of the samples. The ASVs showed significant differences between the intervention groups and the NI group after intervention (adjusted p< 0.05, log2 fold change >2) (Fig. 2c). Compared with the NI group, we found that 15 ASVs were enriched and 5 ASVs decreased in the AED group; and 13 ASVs were enriched in the AEx group, while 9 ASVs were enriched and 6 ASVs were decreased in the Diet group after the intervention. Among these ASVs, ASV2077 (AED: p=0.018, AEx p< 0.01, Diet: p=0.032) and ASV2513 (AED: p =0.012, AEx: p< 0.01, Diet: p=0.013) belonged to Bacteroides, and ASV3942 (AED: p=0.030, AEx: p=0.024, Diet: p< 0.01), which belongs to Ruminococcus, increased in all intervention groups. In addition, ASV5361, which belongs to Lachnospiraceae, increased in both AEx (p=0.027) and AED (p=0.037) groups. ASV2440, belonging to Bacteroides, increased in both Diet (p< 0.01) and AED (p< 0.01) groups. Noticeably, some ASVs from the same family or genus showed different behaviors. For example, ASV 4432 in Lachnospiraceae was enriched in the AED (p< 0.01) group but decreased after the diet (p=0.023) intervention compared with the NI group. To investigate whether changes of microbiome at the ASV level were similar to the changes at the genus level, we further performed a LEfSe analysis of the same parameters on the genus level. We found that those observed changed ASVs, if belonging to the same genus, indeed, they do have a similar trend at the genus level (Fig. S3). We next performed partial Spearman correlation analysis to assess change of ASVs with clinical biomarkers after adjustment for body weight and fat mass (n=46, Fig. 2d). We found that three altered ASVs were significantly negatively associated with the reduction of HFC [ASV2468, belonging to Bacteroides (r=−0.31, p=0.040), ASV3307, belonging to Ruminococcaceae (r=−0.32, p=0.030), and ASV4538, belonging to Lachnospira (r=−0.37, p=0.012)]. In addition, ASV478, which belongs to Phascolarctobacterium (r=−0.32, p=0.033); ASV1715, which belongs to Alistipes (r=−0.35, p=0.022); and ASV 5195 and ASV5305, which belong to Lachnoclostridium (r=−0.30, p=0.045, r=−0.34, p=0.023, respectively), were negatively correlated with HbA1c. Only ASV776, • AED: n=23 • AEx: n=22 • Diet: n=22 • NI: n=18 • Clinical phenotypes (HFC, fat mass, serum biomarkers … ) Baseline to follow-up: 85 pairs of individuals • Gut microbiome: 16S rRNA sequence Baseline to follow-up: 76 pairs of individuals • AED: n=20 • AEx: n=20 • Diet: n=21 • NI: n=15 • Gut microbiome: Metagenomics • SCFA Baseline to follow-up: 42 pairs of individuals • AED: n=12 • AEx: n=12 • Diet: n=10 • NI: n=8 Measurement and Data Collection AED (n=29) AEx (n=29) Diet (n=28) NI (n=29) Randomization (n=115) Baseline AED (n=23) Follow-up AEx (n=22) Diet (n=22) NI (n=18) Intervention Trial (average 8.6M) AED AEx Diet NI 0 10 20 30 40 VO2max (ml/kg/min) 0 20 40 Dietary fiber (g/d) 0M 8.6M * * * * 0M 8.6M 0M 8.6M 0M 8.6M a b 653 2065 143 346 270 1004 57 18 2008 78 2061 519 1034 470 111 336 269 723 1003 178 4014 217 345 6018 4003 4019 224 327 160 22 46 145 14 24 43 3016 26 6002 2064 442 2006 256 2051 592 304 599 1052 70 157 2025 87 354 1010 33 44 417 29 2018 402 2067 661 Metagenomic and 16s sequencing 16s sequencing only 10 -10 AED AEx Diet -20 -30 -40 -50 -60 0 HFC percentage change Non-responderResponder c Fig. 1 Summary of clinical study. a Study design. The measurements were performed before and after intervention, and comparison of dietary fiber intakes and estimated physical fitness indicated by VO2max before and after intervention are given. The white point, box range, line range and density plot width represent median, interquartile range, 95% confidient interval and frequence, respectively. AEx: aerobic exercise; AED: AEx+Diet; NI: no intervention. Comparisons within groups were performed using two-tailed paired t-test. *: Significance of p< 0.05. bResponders and low/non-responders according to change of HFC (the cutoff value was -5%) after intervention. AED: n=20, AEx: n=20, Diet: n=21. cData collection. Body fat mass was measured by DXA (dual X-ray densitometry); hepatic fat content (HFC) was determined by 1HMRS (proton magnetic resonance spectroscopy); fecal samples were collected for microbiome assessments with both 16 s rRNA and metagenome assessments; and blood samples were used for clinical biomarkers assessments. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29968-0 ARTICLE NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications 3 belonging to Erysipelotrichaceae (r=0.36, p=0.016), was positively associated with HBA1c. Although a previous study23 showed the gut microbiota may benefit humans via SCFAs production from carbohydrate fermentation, fecal SCFAs did not change significantly from baseline to follow-up within groups in this experiment. However, we found that some alterations in gut microbes correlated with changes in SCFAs (Fig. 2d). The level of butyric acid significantly correlated with 9 ASVs. Notably, ASV3718 in the genus Faecalibacterium was significantly negatively correlated with 4 out of all 6 SCFAs, while ASV1989 in the genus Bacteroide was positively correlated with 4 SCFAs (Fig. 2d). To assess microbial functions, we also applied shotgun sequencing in a subset of samples (Fig. 1c) and annotated the data to the KEGG pathway (Fig. 3). By using LEfSe analysis (p< 0.05, LDA score > 2), we found that in total, 64 specific pathways were significantly different between the AED/AEx/Diet groups and NI group after interventions, most of which fall into the ‘carbohydrate metabolism’,‘energy metabolism’,‘glycan biosynthesis and metabolism’,‘lipid metabolism’,‘amino acid metabolism’and ‘metabolism of cofactors and vitamins’function pathways. The differences in functional pathways of the three intervention groups relative to the NI group were generally AED AEx Diet NI 0M 8.6M 0M 8.6M 0M 8.6M 0M 8.6M 1 2 3 4 5 Shannon diversity * *** ** ** -0.2 0.0 0.2 -0.5 -0.25 0.00 0.25 PCoA.1 PCoA.2 -4 -2 0 2 4 sign_interv -4 -2 0 2 4 sign_interv -4 -2 0 2 4 sign_interv Adlercreutzia 356 363 Bacteroides 2077 2513 2434 2440 2081 2374 2375 2458 2500 2515 Bifidobacterium 189 Bilophila 82 Collinsella 417 Eggerthella 341 Klebsiella 1415 Lachnospiraceae 4432 5361 Lactobacillus 1025 Megamonas 701 Phascolarctobacterium 478 Prevotella 2284 2306 Pseudomonas 1533 Roseburia 4876 4903 Ruminiclostridium 3378 Ruminococcaceae 3305 3102 3126 Ruminococcus 3942 Sutterella 1585 unclassified 4768 5016 5407 Rich side color AED AEx Diet NI a c b d After intervention Comparision of AED and NI * + + * * + + * * * * * * * * * + * * * * * * * * + + * * * * * * * * * * * * * * * * * + * * * * + * * * * * * * * * * * * + * * * * * * * * * * * * * * * + * * Acetic_acid Propionic_acid Isobutyric_acid Butyric_acid Isovaleric_acid Valeric_acid SCFA AST HBA1c HFC VO2max Energy_kcal Fiber_g Protein_E Fat_E Carbohydrate_E HOMA ISI HDL LDL Trigly [Eubacterium] coprostanoligenes group; ASV_3868 Blautia; ASV_5102 Alistipes; ASV_1764 Alistipes; ASV_1761 UBA1819; ASV_3535 Alistipes; ASV_1715 Phascolarctobacterium; ASV_478 Flavonifractor; ASV_3199 Dialister; ASV_575 Bacteroides; ASV_2515 Collinsella; ASV_417 Lachnospira; ASV_4538 Faecalibacterium; ASV_3579 Lachnospiraceae ND3007 group; ASV_4432 Romboutsia; ASV_1100 Bacteroides; ASV_2480 Subdoligranulum; ASV_3758 unclassify; ASV_4768 Bacteroides; ASV_2468 Ruminococcaceae NK4A214 group; ASV_3307 Faecalibacterium; ASV_3718 Roseburia; ASV_4879 Blautia; ASV_5037 Intestinibacter; ASV_1125 [Eubacterium] hallii group; ASV_4515 Blautia; ASV_5117 Ruminococcaceae UCG 002; ASV_3131 Roseburia; ASV_4847 Erysipelotrichaceae UCG 003; ASV_776 Tyzzerella 3; ASV_4337 Ruminococcus 2; ASV_3436 Bacteroides; ASV_1973 Bacteroides; ASV_1983 Fusicatenibacter; ASV_5165 Anaerostipes; ASV_4591 unclassify; ASV_4886 Lachnoclostridium; ASV_5195 Lachnoclostridium; ASV_5305 Lachnospira; ASV_4523 Lachnoclostridium; ASV_5210 Bacteroides; ASV_2053 [Ruminococcus] gnavus group; ASV_4682 Escherichia Shigella; ASV_1501 Roseburia; ASV_4903 Bacteroides; ASV_1989 Bifidobacterium; ASV_152 Anaerostipes; ASV_4595 ASV Comparision of AEx and NI Comparision of Diet and NI AED AEx Diet NI −0. 4 −0.2 00.2 0.4 Fig. 2 Changes in the structure of gut microbiota before and after intervention and relationship of microbial abundance with clinical parameters. aAlpha diversity (Shannon index) of gut microbiota was significantly decreased in the NI group by contrast to the intervention groups, which maintained their diversity after the intervention. The white point, box range, line range and density plot width represent median, interquartile range, 95% confidient interval and frequence, respectively. ANCOVA for repeated measures (2 factor interactions: group × time) and controlled for change of body weight, baseline value and intervention duration followed by Sidak correction for multiple comparison between the groups. Contrast results (K Matrix) were used to localize the group differences: *p < 0.05, by contrast to the NI group. AED 8.6 m vs NI 8.6 m, p=0.011; AEx 8.6 m vs NI 8.6 m, p=0.007; Diet 8.6 m vs NI 8.6 m, p=0.025; NI 0 m vs NI 8.6 m, p=0.031. AED: n=20, AEx: n=20, Diet: n=21 and NI: n=15. bWeighted Unifrac distance of samples was calculated from the QIIME 2 ASV level. Ellipsoids represent a 50% confidence interval surrounding each group. cDifferential ASVs between the AED/AEx/ Diet groups and NI group (LEfSe analysis, adjusted p< 0.05, log2 fold change > 2). Genus annotations on the left column covers multiple right ASVs. AED: n=20, AEx: n=20, Diet: n=21 and NI: n=15. dHeatmap of the Spearman’s correlation coefficients between change of ASVs as assessed by 16 S and clinical parameters independent of body weight and fat mass. Statistically significant coefficients are marked by * and +, which means p< 0.05 and FDR < 0.1respectively. Only ASVs with significant correlations (at least one based on p value) are shown. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29968-0 4NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications similar, however, the AED group had more differential pathways than other groups. Notably, within ‘lipid metabolism’,‘Sphingolipid metabolism’was more vigorously different from that of other pathways (Fig. 3). In addition, most functional pathways were more abundant in the NI group; only ‘glycan biosynthesis and metabolism’was more abundant in the intervention groups. Taken together, the above findings indicate that exercise or/and diet intervention(s) distinctly alter the abundance and function of microbiome, which are associated with changes in HFC and SCFAs. Intervention induced change in gut microbiota co-occurrence network. It has been shown that bacterial species in the human gut may survive, adapt, and decline as interdependent functional groups (guilds) responding to environmental perturbations24,25. To identify bacteria in the gut ecosystem that responded as functional groups to interventions, we adopted “co-abundance groups (CAGs)”to analyze the community structure in the microbial ecosystem24. We used the SparCC algorithm to calculate the correlation coefficients among 279 ASVs shared by at least 20% of the samples from all the groups and time points26. These 279 ASVs were then clustered into 35 CAGs (Table S2) and their abundance difference between groups were showed in Fig. S4a. Network analysis can disentangle microbial co-occurrence and provide comprehensive insight into the microbial assembly patterns and community structures27. Using Spearman correlation analysis, we constructed a co-occurrence network (Fig. 4a–d) to illustrate the potential interactions among the 35 CAGs in each group with Glycosaminoglycan degradation Proximal tubule bicarbonate reclamation Glycosphingolipid biosynthesis globo and isoglobo series ABC transporters Various types of N glycan biosynthesis Xylene degradation Sulfur metabolism Valine leucine and isoleucine biosynthesis Lipopolysaccharide biosynthesis Pantothenate and CoA biosynthesis Glyoxylate and dicarboxylate metabolism Bacterial secretion system Protein processing in endoplasmic reticulum Neomycin kanamycin and gentamicin biosynthesis Tyrosine metabolism Lysine degradation Sphingolipid metabolism Zeatin biosynthesis Polyketide sugar unit biosynthesis Glycerophospholipid metabolism Lysosome Sulfur relay system Biotin metabolism N Glycan biosynthesis Chloroalkane and chloroalkene degradation Folate biosynthesis Arginine biosynthesis Pyruvate metabolism Methane metabolism Nicotinate and nicotinamide metabolism Limonene and pinene degradation Butanoate metabolism Oxidative phosphorylation Autophagy yeast Glycerolipid metabolism C5 Branched dibasic acid metabolism Chlorocyclohexane and chlorobenzene degradation Tryptophan metabolism beta Alanine metabolism Prion diseases Glycosphingolipid biosynthesis ganglio series Nitrogen metabolism Photosynthesis Phosphotransferase system -PTS Cysteine and methionine metabolism D Arginine and D ornithine metabolism Pentose phosphate pathway Bacterial chemotaxis Purine metabolism Cyanoamino acid metabolism Other glycan degradation Selenocompound metabolism Propanoate metabolism Benzoate degradation Atrazine degradation Fatty acid degradation Naphthalene degradation Proteasome Biosynthesis of unsaturated fatty acids Porphyrin and chlorophyll metabolism Glycolysis / Gluconeogenesis Lipid metabolism Carbohydrate metabolism Glycan biosynthesis and metabolism Amino acid metabolism Energy metabolism Metabolism of cofactors and vitamins Intervention Group AED AEx Diet Fig. 3 Identified pathways by KEGG between intervention groups and NI group after intervention. Only significant pathways with a linear discriminate analysis (LDA) score > 2.0 and p< 0.05 are shown. Six colors indicate six primary function pathways. The bars inside the circle indicate that 3 intervention groups had higher values than the NI group, and the bars outside the circle indicate that the intervention groups had lower values than the NI group. The bar length represents LDA score. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29968-0 ARTICLE NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications 5 baseline and follow-up samples. Each node in the network represents a CAG, and each edge represents a significant correlation (p< 0.05) between two CAGs. We found that the number of edges and average degrees that signify the connectivity of CAGs in the gut microbiota network were higher after intervention than at baseline in the AED, AEx and Diet groups, but not in the NI group (Fig. 4e). In addition, the network robustness28 was significantly decreased in the NI group but maintained or increased in all intervention groups (Fig. 4f). The greatest increase in robustness score was in the AEx group (+29.3%), followed by the Diet (+12.8%) and AED ( +3.1%) groups. Because the complexity caused by the nonlinear dynamics of the nodes is related to the network architecture, we assessed the distribution of the network degree to judge the topological features of the microbial network (Fig. S4b). We found that most of the connections of the co-occurrence network were CAG32 CAG35 CAG25 CAG10 CAG5 CAG21 CAG15 CAG26 CAG8 CAG6 CAG27 CAG28 CAG4 CAG30 CAG13 CAG31 CAG1 CAG34 CAG22 CAG3 CAG18 CAG9 CAG2 CAG7 CAG23 CAG29 CAG17 CAG19 CAG20 CAG16 CAG12 CAG33 CAG24 CAG14 CAG1 CAG15 CAG4 CAG20 CAG16 CAG11 CAG14 CAG7 CAG3 CAG12 CAG10 CAG30 CAG27 CAG22 CAG6 CAG8 CAG35 CAG24 CAG17 CAG29 CAG33 CAG25 CAG5 CAG26 CAG31 CAG19 CAG18 CAG9 CAG28 CAG13 CAG2 CAG21 CAG23 CAG32 CAG17 CAG13 CAG22 CAG31 CAG34 CAG19 CAG1 CAG14 CAG20 CAG3 CAG24 CAG11 CAG23 CAG7 CAG5 CAG16 CAG9 CAG32 CAG10 CAG15 CAG27 CAG33 CAG4 CAG35 CAG28 CAG30 CAG29 CAG21 CAG18 CAG12 CAG2 CAG6 CAG26 CAG8 CAG25 CAG7 CAG2 CAG31 CAG1 CAG18 CAG29 CAG8 CAG27 CAG3 CAG22 CAG30 CAG21 CAG12 CAG6 CAG25 CAG23 CAG5 CAG17 CAG20 CAG11 CAG19 CAG24 CAG9 CAG33 CAG28 CAG13 CAG16 CAG14 CAG26 CAG34 CAG15 CAG10 CAG35 CAG32 CAG4 a c b d 0 50 100 Edge Number 02 4 Mean Degree AED AEx Diet NI 56 58 52 71 83 96 69 54 3.294 4.000 3.714 4.176 4.882 5.486 3.943 3.600 e Remove nodes percentage AED AEx Diet NI 0.0 0.5 1.0 Robustness score=11.167 score=11.517 score=10.951 score=14.162 score=12.479 score=14.081 score=12.963 score=10.925 0.5 100.5 100.5 10 0.5 10 f Fig. 4 Co-occurrence networks before and after intervention in paired participants. a–dCo-occurrence networks of CAGs (co-abundance groups) before and after intervention in different groups (a: AED: n=20; b: AEx: n=20; cDiet : n=21; d: NI: n=15). The nodes represent CAGs and the edges represent statistically significant differences (p< 0.05) via Spearman correlation between CAGs without adjustment. The style of the edge is set to a dashed line and a solid line, which represent negative and positive correlation, respectively. Gray edges indicate correlation before intervention; colored edges indicate correlation after intervention. eNetwork connection properties. fNetwork robustness test. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29968-0 6NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications concentrated in a few nodes, which is an important features of a scale-free network. This kind of network with a power-law distribution comprise highly connected nodes, which are defined as hubs29. The hubs in a microbial network have been proposed as keystone taxa, as their removal has been computationally shown to cause a drastic shift in the composition and functioning of a microbiome20. In our data, the nodes (CAG6, CAG7, CAG8, and CAG28) had more than 4% of the total number of connections in each network. Therefore, these CAGs can be regarded as hubs and maybe the keystone taxa of the gut microbial community of each group (Fig. S4c). Interestingly, the abundance of these four CAGs increased significantly after intervention in the AED group compared with the NI group, but did not change significantly in the AEx and Diet groups (Fig. S4a). The most abundant ASVs among these four CAGs belong to Bacteroides,Ruminococcaceae,Alistipes and Subdoligranulum. Collectively, these findings suggest that interventions improve the stability of the microbial ecosystem by reconstructing the network among keystone taxa. Personalized gut microbial networks can predict the intervention efficiency of individual patients. Since not all subjects responded to the intervention in a similar manner (Fig. 1b), we assessed whether these response variations were associated with the individual characteristics of the gut microbiota. First, we used differentially altered ASVs abundance in each intervention group (Fig. 2c) to classify the responders and low/non-responders, we found that only the baseline ASVs abundance in the Diet group showed the prediction power with a receiver operating characteristic (ROC) area under the curve (AUC) of 0.65 (specificity =0.58, sensitivity =0.67). No significant prediction power for the classification in the AED group (AUC =0.53, specificity =0.40, sensitivity =0.67) and the AEx group (AUC =0.52, specificity =0.50, sensitivity =0.60) was found. Second, since ASVs abundance could not describe the individual network characteristics, we developed a Single SparCC network method by combining a SparCC network26 with a sample specific network (SSN)30 (Fig. 5a). We used baseline metagenomics data to build a network for each patient (Fig. 5b) and found that responders in all three intervention groups tended to have more interactions between species than low/non-responders, but no statistical significance was observed (Fig. S5). Furthermore, in the linear regression analysis, the edge number in the AEx group significantly predicted the change of HFC (Fig. 5c). Then we performed an unsupervised classifier to differentiate the responders from the low/non-responders by each Single SparCC network edge number. We found that AUC was about 70% in different intervention groups (Fig. 5d), and the AUC of supervised Least absolute shrinkage and selection operator (LASSO) classifiers were higher in all three groups (Fig. 5e). For the purpose of comparison, we also tested the ability of age, body weight (WT) and body mass index (BMI) to distinguish the responsiveness of the interventions (Fig. S6). The result showed that these clinical parameters were not able to differentiate the responders from the low/non-responders (age: AUC =0.437, p=0.347; WT: AUC =0.470, p=0.655; BMI: AUC =0.519, p=0.771), except for BMI in the AED group. Considering the small sample size of metagenomic data, we performed the same Single SparCC network analysis on the intervention groups’samples by using the 16 S rRNA gene sequencing data (Fig. S7). Consistent with the metagenomic results, the responders of the three intervention groups tended to have more network edges than low/non-responders, but none of them were statistically significant (Fig. S8a). In regression analysis, we found that, in addition to the AEx group, the baseline edge numbers in the AED group also showed significant correlation with the HFC change; however, this was not observed in the Diet group (Fig. S8b). For supervised LASSO classifiers, ROC showed that the edge number of ASV data used to predict the responder exceeded 70% AUC in all three intervention groups (Fig. S8c), which was higher than the metagenomics data. However, the effect of unsupervised classifiers deteriorated and was almost ineffective in the Diet group (Fig. S1d). These results indicated that our Single SparCC networks could be used to construct personalized gut microbial networks for each individual sample and thus differentiate the responders from the low/nonresponders, particularly for the exercise intervention. Taken together, these analyses demonstrate that the individual baseline gut microbial network can predict the response to exercise intervention. Discussion In the current study, in patients with NAFLD and pre-diabetes, we identified the characteristics of the gut microbiota responding to an 8.6-month aerobic exercise and/or low-carbohydrate dietary intervention. We showed that after combined exercise and diet intervention, changes are prominent in gut microbial composition in which keystone taxa became divergent but their connections remained stable. By contrast, with exercise or diet intervention alone, effects were significant regarding network connectivity and robustness between taxa. Moreover, the personalized microbial network was able to predict the intervention efficiency of individual patients. Members from Ruminococcus have been reported to produce SCFAs from complex carbohydrates31,32. A recent study33 showed that Ruminococcaceae was negatively correlated with the fibrosis severity. In agreement with this, members in Ruminococcaceae was found negatively correlation with HFC in our study. In addition, some members of Bacteroides contribute to the release of energy from dietary fiber and starch, and they are likely to be a major source of propionate34. Accordingly, we found that members in Bacteroides was positively correlated with propionic acid, isobutyric acid and isovaleric acid. A previous study showed that compared to western countries, the abundance of Bacteroides in NAFLD is much lower in Chinese individuals35. However, our result suggested an important energy-extracting role of members in Bacteroides, although their abundance is relatively low in Chinese population. Of note, there are currently conflicting results about the shifts of microbiota at the taxon level in NAFLD-related studies. For example, some studies showed that Bacteroides abundance was higher in NASH (nonalcoholic steatohepatitis, a stage of NAFLD progression) than in non-NASH36. Another article showed the phylum Bacteroidetes (44.63%) tended to be more abundant in healthy subjects than in NAFLD35. Our results showed that members in Bacteroides abundance was increased after both diet and exercise intervention, and this increase was correlated with decreased HFC. However, previous studies showed higher fiber intake was marginally associated with lower abundance of Bacteroides uniformis37,andBacteroides was increased in participants with obesity but decreased in lean participants after exercise intervention38. These inconclusive results may be ascribed to differences in ethnicity, living environment and lifestyle among the cohorts. Thus, the mechanism of action of the key microbiota needs to be investigated in future studies. The long-term adaptation after intervention may account for maintenance of gut microbiota alpha diversity in all intervention groups in contrast to the NI group in which alpha diversity decreased. However, the connectivity and robustness of the microbiota co-occurrence networks improved in three NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29968-0 ARTICLE NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications 7 intervention groups. Various ecological microbiome studies39,40 have shown that healthy people have a higher average network degree and hence greater connectivity of their gut microbiota network41. A poorly developed microbial network usually has lower functionality due to fewer taxa present that can support their function in ecosystems39. In combination with these findings, our results indicate that microbiota as an ecosystem were more healthy and stable after intervention. Interestingly, when Responders Non-responders AED AEx Diet R2=0.26 P = 0.0914 R2=0.45 P = 0.0166 R2=0.35 P =0.0739 -60 -40 -20 0 3000 9000 edge number HFC change Group AED AEx Diet Specificity Sensitivity 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 AUC AED 0.90 AEx 0.94 Diet 0.76 6000 120000 b cd 10 1 1. Importing species abundance of N samples 2N-1 N Species A 08… 100 0 0 Species B 10… 1000 Species C 100 2…830 3. Distinguishing responders, non/low-responders by Single SparCC network Responders Non/low-responders 2. Estimating Single SparCC network for each sample Species A Species B Estimating edges in Nth samples by comparing their abundance to the N-1 distribution Estimating the distribution of species abundance in N-1 samples Single SparCC network for each sample e Specificity Sensitivity 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 AUC AED 0.77 AEx 0.74 Diet 0.72 a Fig. 5 Intervention efficiency predicted by individual baseline gut microbial networks on metagenomics species data. a Prediction of responders using our single SparCC network method. bThe gut microbial network for each individual before intervention based on metagenomics species data by the Single SparCC network method. Three numbers under every network represent subject ID, HFC change and edge numbers. cLinear regressions of edge number (Single SparCC network) with the change of HFC after intervention. AED: n=12, AEx: n=12, Diet: n=10. dROC curve of edge number for differential responders from low/non-responders within each intervention group by unsupervised classification. AED: n=12, AEx: n=12, Diet: n=10. eROC curve by supervised Lasso model performance. AED: n=12, AEx: n=12, Diet: n=10. ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29968-0 8NATURE COMMUNICATIONS | (2022) 13:2555 | https://doi.org/10.1038/s41467-022-29968-0 | www.nature.com/naturecommunications