ARTICLE

[Link] OPEN

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 ✉,
1234567890():,;

Luonan Chen4,11,12 ✉ & Sulin Cheng 1,2,5 ✉

Exercise and diet are treatments for nonalcoholic fatty liver disease (NAFLD) and pre-
diabetes, 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 robust-
ness 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.

A full list of author affiliations appears at the end of the paper.

NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications 1

ARTICLE NATURE COMMUNICATIONS | [Link]

N
onalcoholic fatty liver disease (NAFLD) is the most Results
common chronic liver disease with prevalence estimates Participant characteristics and gut microbiota research design.
ranging from 25% to 45% in worldwide1,2. NAFLD is This study was an 8.6-month, four-arm, randomized trial (Fig. 1a,
closely associated with type 2 diabetes (T2D)3, most likely via a Fig. S1) and participant characteristics have been reported in our
common pathophysiological mechanism - insulin resistance4. previous publication21. Briefly, 115 participants were recruited from
Therefore, interventions targeting the synergistically pathogenic 7 health clinical service centers in the Shanghai Yangpu district.
mediator of NAFLD and diabetes are likely to be a rational They were randomized into four groups: aerobic exercise inter-
approach in the prevention and treatment of the comorbidity vention (AEX, n = 29), fiber-enriched low-carbohydrate diet
condition. intervention (Diet, n = 28), aerobic exercise combined with diet
The gut microbiota, which contributes to the metabolic health intervention (AED, n = 29) and no intervention without guided
of the human host, may work as one of the targets for the exercise and dietary intake (NI, n = 29). Of all participants, 85
treatment of NAFLD. Emerging data demonstrate that dysbiosis individuals completed the intervention trial.
of the gut microbiota is linked with NAFLD and T2D5,6. For the primary outcome of liver fat content, we previously
Although the exact mechanism(s) require further elucidation, found that hepatic fat content (HFC) was significantly reduced in
inflammation, damage to the intestinal membrane, and translo- the exercise AEx (–24.4%), diet (–23.2%), and AED (–47.9%)
cation of bacteria have all been suggested7. Drugs with antifibrotic groups by contrast to the 20.9% increase in the NI group (p < 0.001
or anti-inflammatory treatments for various stages of NAFLD in for all) after intervention21. Of note, 91% of the subjects in the
human trials is pending8. Currently, increased physical activity AED group decreased their HFC, and the corresponding figures
and dietary modifications are the only effective therapeutic were 68% in the AEx group, and 86% in the Diet group. In
options for NAFLD management, and the mechanisms of these contrast, 72% in the NI group increased their HFC during the
interventions have been associated with modulation of gut intervention period. However, for the primary outcome of glucose
microbiota and its metabolites9. For example, a low-carbohydrate metabolism, no significant remission or progression of prediabetes
diet (LCD) has been reported to improve fatty liver metabolism was found between the intervention and NI groups based on the
and to promote rapid shifts of the gut microbiota composition in glycated hemoglobin A1c (HbA1c)21. These results indicated that
NAFLD patients10. Regular aerobic exercise has been shown to our intervention was effective mainly for HFC reduction.
reduce hepatic fat in obesity11,12, and some studies have shown We found that not all the subjects responded to the intervention
that exercise increases gut microbial diversity13 and alters the in a similar manner. Therefore, we stratified the participants
composition and functional capacity of gut microbiota14,15. into responders (HFC decreased more than 5%) and low/non-
However, the underlying mechanism of the beneficial effects of responders (HFC decrease less than 5% or increased) (Fig. 1b). For
exercise and diet on NAFLD, as well as the related metabolic the primary outcome of gut microbiota composition, seventy-six
disorders through regulation of the gut microbiota, still remain to subjects provided paired stool samples at baseline and after the
be elucidated. intervention for 16 s rRNA gene sequencing. The demographics and
One of the challenges in constructing the relationship clinical variables of those who have had gut microbiota results are
between the gut microbiota and improvement of NAFLD during presented to Table S1. In addition, to better understand the
exercise and/or dietary interventions is that not all patients intervention impact on the function of microbiome, we selected a
respond to the intervention in a similar manner. Earlier studies subset of the cohort (n = 42) with the best or worst response, in
have shown that there was an individual variation in response to terms of HFC reduction, to the intervention, and analyzed gut
exercise, with some subjects experiencing greater improvement microbiota by the metagenomics data (Fig. 1c).
than others16. The low/non-responders could be up to 50% in The microbiota data and clinical parameters were used to (1)
terms of change of gut microbiota after exercise intervention17. characterize the change of gut microbiota composition in
Some researchers have suggested that the variation in respon- response to interventions; (2) determine associations of the
siveness might be due to genotypic and phenotypic factors18,19. microbiome with the level of physical fitness, fat mass, serum
Thus, large inter-individual variance may mask the changes of biomarkers and short-chain fatty acids (SCFAs); (3) analyze the
the microbiome in studies using traditional statistical analytical metabolic function shift of the microbiome in response to
methods. Currently, it is widely accepted that, in the microbial interventions; (4) discover the co-occurrence network features of
community, the keystone taxa are drivers of microbiome the microbiome ecosystem before and after intervention; and (5)
structure and function, and in particular, their interaction net- predict personalized response to intervention based on the
work, which plays an important role in microbial functions and baseline gut microbiota composition and network.
disease progression20. Therefore, to establish an effective inter-
vention strategy, it may be worthwhile to analyze the microbiota
correlation network at the microbial community level, while Changes of structure and function of gut microbiota induced
assessing the physiological mechanisms at the individual level by exercise and/or dietary intervention of participants with
and further exploring individual microbial networks that comorbidity of NAFLD and prediabetes. We first performed
underlie the differences between responders and low/non- 16 S rRNA gene sequencing of fecal samples collected before and
responders of various interventions. after interventions and obtained a total of 5421 amplicon
In the present study, we analyzed the composition and meta- sequence variants (ASVs). We found that the alpha diversity
bolic pathways of the gut microbiota, constructed a co-occurrence (Shannon index) of gut microbiota was significantly decreased in
network at the population level, and developed personalized gut the NI group (p = 0.045) with large individual variance by con-
microbial networks. This individual network enabled us to trast to the intervention groups, which maintained their diversity
identify the microbial signature and interaction of taxa at a per- after the intervention (Fig. 2a, NI vs. AED p = 0.011, NI vs. AEx
sonal resolution, and to further differentiate responders from low/ p = 0.007 and NI vs. Diet p = 0.025, respectively, analysis of
non-responders after intervention. Taken together, our study with variance with repeated measures adjusted for change of body
relatively long-term intervention provides both statistical and weight, baseline value and intervention duration). The changes in
sample-specific network insights into the complex gut microbial microbial diversity may reflect the natural progression of NAFLD
ecosystem in patients with NAFLD and glucose metabolism in the NI group. Moreover, the principal coordinates analysis
impairment. (PCoA) based on weighted UniFrac distances shows the

2 NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications

NATURE COMMUNICATIONS | [Link] ARTICLE

a
AED AEx Diet NI
Intervention Trial (average 8.6M) 40 * *

VO2max (ml/kg/min)
Randomization (n=115) 30

20
Baseline

AED (n=29) AEx (n=29) Diet (n=28) NI (n=29) 10

0
40 * *

Dietary fiber (g/d)

20
Follow-up

AED (n=23) AEx (n=22) Diet (n=22) NI (n=18)

0

0M 8.6M 0M 8.6M 0M 8.6M 0M 8.6M

Metagenomic and 16s sequencing c

b 16s sequencing only Measurement and Data Collection

Non-responder
AED AEx Diet
10 • Clinical phenotypes • AED: n=23
(HFC, fat mass, serum biomarkers … ) • AEx: n=22
0 Baseline to follow-up: • Diet: n=22
HFC percentage change

85 pairs of individuals • NI: n=18

-10
• Gut microbiome: 16S rRNA sequence • AED: n=20
-20 Baseline to follow-up: • AEx: n=20
76 pairs of individuals • Diet: n=21
-30

Responder
• NI: n=15
-40 • Gut microbiome: Metagenomics • AED: n=12
• SCFA • AEx: n=12
-50 Baseline to follow-up: • Diet: n=10
42 pairs of individuals • NI: n=8
-60
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

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. b Responders 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. c Data 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.

significant difference between NI and the other groups after the increased in all intervention groups. In addition, ASV5361,
interventions (p < 0.05, PERMANOVA, Fig. 2b). These results which belongs to Lachnospiraceae, increased in both AEx
indicate that microbial diversity deteriorates with the increased (p = 0.027) and AED (p = 0.037) groups. ASV2440, belonging
HFC, while exercise and dieting may help to maintain the to Bacteroides, increased in both Diet (p < 0.01) and AED
diversity of the gut microbiota. To avoid only included the rela- (p < 0.01) groups. Noticeably, some ASVs from the same family
tive abundance data in the analysis which may introduce the bias, or genus showed different behaviors. For example, ASV 4432 in
we then measured the absolute abundance of microbiome by Lachnospiraceae was enriched in the AED (p < 0.01) group but
quantitative real-time PCR (qPCR) targeting the 16 S rRNA gene. decreased after the diet (p = 0.023) intervention compared with
We found that there were no significant differences between the NI group. To investigate whether changes of microbiome at
baseline and after intervention as well as among the groups the ASV level were similar to the changes at the genus level, we
(Fig. S2a). Moreover, the total bacterial content of each subject in further performed a LEfSe analysis of the same parameters on the
the same group at two time points also did not differ significantly genus level. We found that those observed changed ASVs, if
(Fig. S2b–e). belonging to the same genus, indeed, they do have a similar trend
Furthermore, we performed Linear discriminant analysis Effect at the genus level (Fig. S3).
Size (LEfSe) analysis22 of ASVs that appeared in more than 20% We next performed partial Spearman correlation analysis to assess
of the samples. The ASVs showed significant differences between change of ASVs with clinical biomarkers after adjustment for body
the intervention groups and the NI group after intervention weight and fat mass (n = 46, Fig. 2d). We found that three altered
(adjusted p < 0.05, log2 fold change >2) (Fig. 2c). Compared with ASVs were significantly negatively associated with the reduction of
the NI group, we found that 15 ASVs were enriched and 5 ASVs HFC [ASV2468, belonging to Bacteroides (r = −0.31, p = 0.040),
decreased in the AED group; and 13 ASVs were enriched in ASV3307, belonging to Ruminococcaceae (r = −0.32, p = 0.030),
the AEx group, while 9 ASVs were enriched and 6 ASVs were and ASV4538, belonging to Lachnospira (r = −0.37, p = 0.012)].
decreased in the Diet group after the intervention. Among these In addition, ASV478, which belongs to Phascolarctobacterium
ASVs, ASV2077 (AED: p = 0.018, AEx p < 0.01, Diet: p = 0.032) (r = −0.32, p = 0.033); ASV1715, which belongs to Alistipes
and ASV2513 (AED: p = 0.012, AEx: p < 0.01, Diet: p = 0.013) (r = −0.35, p = 0.022); and ASV 5195 and ASV5305, which belong
belonged to Bacteroides, and ASV3942 (AED: p = 0.030, AEx: to Lachnoclostridium (r = −0.30, p = 0.045, r = −0.34, p = 0.023,
p = 0.024, Diet: p < 0.01), which belongs to Ruminococcus, respectively), were negatively correlated with HbA1c. Only ASV776,

NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications 3

ARTICLE NATURE COMMUNICATIONS | [Link]

a b
*
** **
* AED

* AEx
Diet
0.2
5 NI

PCoA.2
Shannon diversity

0.0
3

-0.2
1

AED AEx Diet NI

-0.5 -0.25 0.00 0.25
0M 8.6M 0M 8.6M 0M 8.6M 0M 8.6M PCoA.1

c d
After intervention
* [Eubacterium] coprostanoligenes group; ASV_3868
ASV Comparision of AED and NI Comparision of AEx and NI Comparision of Diet and NI * * Blautia; ASV_5102
Adlercreutzia 356 *Alistipes; ASV_1764
363 *Alistipes; ASV_1761
Bacteroides 2077 * UBA1819; ASV_3535
2513 * * * * Alistipes; ASV_1715
* * * Phascolarctobacterium; ASV_478
2434
* * Flavonifractor; ASV_3199
2440
* Dialister; ASV_575
2081 * Bacteroides; ASV_2515
2374 + Collinsella; ASV_417
2375 * Lachnospira; ASV_4538
2458 * Faecalibacterium; ASV_3579
2500 * Lachnospiraceae ND3007 group; ASV_4432
2515 * * Romboutsia; ASV_1100
Bifidobacterium 189 * * Bacteroides; ASV_2480
Bilophila 82 + Subdoligranulum; ASV_3758
Collinsella 417 * * * unclassify; ASV_4768
Eggerthella 341 * * Bacteroides; ASV_2468
Klebsiella 1415 * * Ruminococcaceae NK4A214 group; ASV_3307
Lachnospiraceae 4432 + + * + + Faecalibacterium; ASV_3718
5361 * * Roseburia; ASV_4879
+ + + Blautia; ASV_5037
Lactobacillus 1025
* * Intestinibacter; ASV_1125
Megamonas 701
* [Eubacterium] hallii group; ASV_4515
Phascolarctobacterium 478
* Blautia; ASV_5117
Prevotella 2284 * * Ruminococcaceae UCG 002; ASV_3131
2306 * Roseburia; ASV_4847
Pseudomonas 1533 * * * Erysipelotrichaceae UCG 003; ASV_776
Roseburia 4876 * Tyzzerella 3; ASV_4337
4903 * * Ruminococcus 2; ASV_3436
Ruminiclostridium 3378 * Bacteroides; ASV_1973
Ruminococcaceae 3305 * Bacteroides; ASV_1983
3102 * * Fusicatenibacter; ASV_5165
3126 * Anaerostipes; ASV_4591
Ruminococcus 3942 + + unclassify; ASV_4886
Sutterella 1585 * Lachnoclostridium; ASV_5195
unclassified 4768 * * Lachnoclostridium; ASV_5305
5016 * * Lachnospira; ASV_4523
5407 * Lachnoclostridium; ASV_5210
* Bacteroides; ASV_2053
-4 -2 0 2 4 -4 -2 0 2 4 -4 -2 0 2 4 * [Ruminococcus] gnavus group; ASV_4682
* * Escherichia Shigella; ASV_1501
sign_interv sign_interv sign_interv * * Roseburia; ASV_4903
* * * * Bacteroides; ASV_1989
* Bifidobacterium; ASV_152
Rich side color * Anaerostipes; ASV_4595
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

AED
AEx
Diet
−0.4 −0.2 0 0.2 0.4
NI

Fig. 2 Changes in the structure of gut microbiota before and after intervention and relationship of microbial abundance with clinical parameters.
a Alpha 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. b Weighted Unifrac distance of samples was
calculated from the QIIME 2 ASV level. Ellipsoids represent a 50% confidence interval surrounding each group. c Differential 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. d Heatmap 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.

belonging to Erysipelotrichaceae (r = 0.36, p = 0.016), was positively To assess microbial functions, we also applied shotgun
associated with HBA1c. Although a previous study23 showed the gut sequencing in a subset of samples (Fig. 1c) and annotated the
microbiota may benefit humans via SCFAs production from data to the KEGG pathway (Fig. 3). By using LEfSe analysis
carbohydrate fermentation, fecal SCFAs did not change significantly (p < 0.05, LDA score > 2), we found that in total, 64 specific
from baseline to follow-up within groups in this experiment. pathways were significantly different between the AED/AEx/Diet
However, we found that some alterations in gut microbes correlated groups and NI group after interventions, most of which fall into
with changes in SCFAs (Fig. 2d). The level of butyric acid the ‘carbohydrate metabolism’, ‘energy metabolism’, ‘glycan
significantly correlated with 9 ASVs. Notably, ASV3718 in the biosynthesis and metabolism’, ‘lipid metabolism’, ‘amino acid
genus Faecalibacterium was significantly negatively correlated with 4 metabolism’ and ‘metabolism of cofactors and vitamins’ function
out of all 6 SCFAs, while ASV1989 in the genus Bacteroide was pathways. The differences in functional pathways of the three
positively correlated with 4 SCFAs (Fig. 2d). intervention groups relative to the NI group were generally

4 NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications

NATURE COMMUNICATIONS | [Link] ARTICLE

Sulfur relay system

Limonen

me
Polyk

endoplasmic reticulum

TS
Zea

Proteaso

te m
m -P
Au

etide
tin b

e rs
e and pin

sys
top
Ly

or t
syste
iosy
hag

su g a

sm
so

p
io

s
so
Ne

oli
nthe

ran
yy

cret
me

i sm
rase
r unit
om

eta
ene deg
eas

Ct
l se
sis
yc

l
Ba

o
dm
AB

ab
ansfe
in

ism
t
ct

biosy

teria
er

et
ka
Pr

un

ol
m
ia
io

Protein processing in
n

ab
lc

o
radation
n m

am

photr

e
Bac
lis

et
di
he

hin
n

om
bo

m
se

yc

th
m
Pr

nit
as a

id
et

in

esis

oc
ot
ox

Phos

ac
es m

or
im
ax
m lis

an

len
al e bo

o
is

D
d
in

in
tub ta

Se
Pu n

d
ge

m
ule la me

an
rin sis

oa
A

nt
em r
bic ta lfu the

an
ine
am
eta ar be Su yn

Cy
s

gin
bo bo

ici
lism na oto
Ph

Ar
te tion

bio
re yla

D
cla hor

sy
ma sp
pho

nt
tio tive

he
n ida olism
etab

sis
Buta Ox en m
g
noate
meta Nitro
sm
bolis Carbohydrate metabolism taboli
m ne me
C5 Bran Glycan biosynthesis and metabolism Metha
ched dib
asic acid ism
metabolis Amino acid metabolism id metabol
m Sphingolip
Glycolysis / Gluconeoge Metabolism of cofactors and vitamins
nesis
Lipid metabolism Glycerophospholipid metabolism
ylat e metabolism
Glyoxylate and dicarbox Energy metabolism Glycerolipid me
tabolism
p a th w ay
te
phospha Bios Fatty a
Pentose s ynth
esis
c id deg
radatio
m rie Por of un n
bolis se p satu
a te meta e s b o hyr
in a rated
ano eri oglo Pa n fa
Prop os
Ni
co nto d chlo tty a
lism
Ch

ion gli nd i
s the rop cids
t abo dat n tin na hyll
lor

m e ra g a a ate te
ate deg is o an an m e
ocy hloro

r u v n e s b dC tab
Py
gly
ca nth glo dn
ico oA
olis
m
clo

ino sy is
C

es tin bio
am bio syn
sis
he

o s d t h am
c pi n
is

the
Gly oli sy he ide
xan
hes

sis
sis

ing bio nt me
alka

h y Fo
ea
the

p pid os lat tab

ynt

os oli bi oli
yc e
ne
nd
syn

olism

Gl ing e bio sm
ios

h r id
chl

sy
and
is

bio

p a Bi
Na
eb
on

os ch nt
es

ot
oro
metab

yc he
Xy

ac
ati

ph

in
th

an

Gl
cin

chlo

s sis
yn

len

m
ad

tha
be

y
lyc

l et
leu
os

po
gr

ab
nze

len
roa
Ng

po
bi

de
de

iso

ionine

Li ol
Benz
an

ed

gr

is
ne
lken
an

esis
of

Tryptophan metabolism

m
lyc

ad
and

Atrazin

eg
Tyrosine metab
lyc

de
es

at
G

ed
oate

rad
d meth
ynth

i
rg

gra
ation

on
N

yp

ne

egr
he

a
st

uci

tio
bios

e degra

degr
Ot

atio
ou

ada
ad

n
e le

ine an
ri

n
Lysine degr
Va

ine

adat

tion
in
Val

Argin

dation
Cyste

ion
oli

Intervention Group
sm

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.

similar, however, the AED group had more differential pathways groups (guilds) responding to environmental perturbations24,25.
than other groups. Notably, within ‘lipid metabolism’, ‘Sphingo- To identify bacteria in the gut ecosystem that responded as
lipid metabolism’ was more vigorously different from that of functional groups to interventions, we adopted “co-abundance
other pathways (Fig. 3). In addition, most functional pathways groups (CAGs)” to analyze the community structure in the
were more abundant in the NI group; only ‘glycan biosynthesis microbial ecosystem24. We used the SparCC algorithm to calculate
and metabolism’ was more abundant in the intervention groups. the correlation coefficients among 279 ASVs shared by at least
Taken together, the above findings indicate that exercise or/and 20% of the samples from all the groups and time points26. These
diet intervention(s) distinctly alter the abundance and function of 279 ASVs were then clustered into 35 CAGs (Table S2) and their
microbiome, which are associated with changes in HFC and abundance difference between groups were showed in Fig. S4a.
SCFAs. Network analysis can disentangle microbial co-occurrence and
provide comprehensive insight into the microbial assembly patterns
Intervention induced change in gut microbiota co-occurrence and community structures27. Using Spearman correlation analysis,
network. It has been shown that bacterial species in the human we constructed a co-occurrence network (Fig. 4a–d) to illustrate the
gut may survive, adapt, and decline as interdependent functional potential interactions among the 35 CAGs in each group with

NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications 5

ARTICLE NATURE COMMUNICATIONS | [Link]

a CAG17 CAG9 b CAG13

CAG29
CAG11
CAG10
CAG16
CAG33 CAG2 CAG3 CAG15

CAG20 CAG6
CAG21 CAG14 CAG12 CAG16
CAG12 CAG35
CAG34 CAG4
CAG25 CAG7
CAG33
CAG26 CAG28 CAG9
CAG1 CAG17 CAG18
CAG5 CAG23
CAG5 CAG20 CAG30
CAG4
CAG27 CAG14 CAG2
CAG24 CAG3 CAG8 CAG27 CAG1 CAG35
CAG22

CAG6 CAG8 CAG24 CAG26

CAG31 CAG25
CAG18 CAG22
CAG21
CAG30
CAG10 CAG7 CAG29
CAG28 CAG31
CAG32 CAG23

CAG19
CAG15 CAG19
CAG13 CAG32

c CAG9 CAG23 d CAG35

CAG12 CAG15
CAG26

CAG29
CAG17 CAG3 CAG18 CAG4
CAG26 CAG19 CAG17
CAG13 CAG8

CAG34 CAG12 CAG7

CAG29 CAG8
CAG20 CAG16
CAG25 CAG27
CAG28 CAG7 CAG23 CAG5

CAG33 CAG25 CAG28

CAG33
CAG14 CAG14 CAG22 CAG3
CAG24 CAG20
CAG19
CAG16 CAG22
CAG24
CAG4 CAG6 CAG11
CAG35 CAG31 CAG27
CAG32 CAG11
CAG30 CAG21 CAG18
CAG1
CAG15 CAG30
CAG6 CAG9
CAG10 CAG2 CAG13
CAG34
CAG2
CAG32
CAG5
CAG31 CAG10 CAG1
CAG21

e AED 56 3.294 f
AED AEx Diet NI
58 4.000
AEx 52 3.714 1.0 score=11.167 score=10.951 score=12.479 score=12.963
71 4.176 score=11.517 score=14.162 score=14.081 score=10.925
Robustness

Diet 83 4.882
96 5.486 0.5
NI 69 3.943
54 3.600
0.0
0 50 100 0 2 4
Edge Number Mean Degree 0 0.5 1 0 0.5 1 0 0.5 1 0 0.5 1
Remove nodes percentage

Fig. 4 Co-occurrence networks before and after intervention in paired participants. a–d Co-occurrence networks of CAGs (co-abundance groups) before
and after intervention in different groups (a: AED: n = 20; b: AEx: n = 20; c Diet : 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. e Network connection properties. f Network robustness test.

baseline and follow-up samples. Each node in the network (Fig. 4f). The greatest increase in robustness score was in the AEx
represents a CAG, and each edge represents a significant correlation group (+29.3%), followed by the Diet (+12.8%) and AED ( + 3.1%)
(p < 0.05) between two CAGs. We found that the number of edges groups.
and average degrees that signify the connectivity of CAGs in the gut Because the complexity caused by the nonlinear dynamics of
microbiota network were higher after intervention than at baseline the nodes is related to the network architecture, we assessed the
in the AED, AEx and Diet groups, but not in the NI group (Fig. 4e). distribution of the network degree to judge the topological
In addition, the network robustness28 was significantly decreased in features of the microbial network (Fig. S4b). We found that
the NI group but maintained or increased in all intervention groups most of the connections of the co-occurrence network were

6 NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications

NATURE COMMUNICATIONS | [Link] ARTICLE

concentrated in a few nodes, which is an important features of a baseline edge numbers in the AED group also showed significant
scale-free network. This kind of network with a power-law correlation with the HFC change; however, this was not observed
distribution comprise highly connected nodes, which are defined in the Diet group (Fig. S8b). For supervised LASSO classifiers,
as hubs29. The hubs in a microbial network have been proposed ROC showed that the edge number of ASV data used to predict
as keystone taxa, as their removal has been computationally the responder exceeded 70% AUC in all three intervention groups
shown to cause a drastic shift in the composition and functioning (Fig. S8c), which was higher than the metagenomics data.
of a microbiome20. In our data, the nodes (CAG6, CAG7, CAG8, However, the effect of unsupervised classifiers deteriorated and
and CAG28) had more than 4% of the total number of was almost ineffective in the Diet group (Fig. S1d). These results
connections in each network. Therefore, these CAGs can be indicated that our Single SparCC networks could be used to
regarded as hubs and maybe the keystone taxa of the gut construct personalized gut microbial networks for each individual
microbial community of each group (Fig. S4c). Interestingly, the sample and thus differentiate the responders from the low/non-
abundance of these four CAGs increased significantly after responders, particularly for the exercise intervention. Taken
intervention in the AED group compared with the NI group, but together, these analyses demonstrate that the individual baseline
did not change significantly in the AEx and Diet groups (Fig. S4a). gut microbial network can predict the response to exercise
The most abundant ASVs among these four CAGs belong to intervention.
Bacteroides, Ruminococcaceae, Alistipes and Subdoligranulum.
Collectively, these findings suggest that interventions improve the
stability of the microbial ecosystem by reconstructing the network Discussion
among keystone taxa. 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
Personalized gut microbial networks can predict the inter- intervention. We showed that after combined exercise and diet
vention efficiency of individual patients. Since not all subjects intervention, changes are prominent in gut microbial composi-
responded to the intervention in a similar manner (Fig. 1b), we tion in which keystone taxa became divergent but their connec-
assessed whether these response variations were associated with tions remained stable. By contrast, with exercise or diet
the individual characteristics of the gut microbiota. First, we used intervention alone, effects were significant regarding network
differentially altered ASVs abundance in each intervention group connectivity and robustness between taxa. Moreover, the perso-
(Fig. 2c) to classify the responders and low/non-responders, we nalized microbial network was able to predict the intervention
found that only the baseline ASVs abundance in the Diet group efficiency of individual patients.
showed the prediction power with a receiver operating char- Members from Ruminococcus have been reported to produce
acteristic (ROC) area under the curve (AUC) of 0.65 (specifi- SCFAs from complex carbohydrates31,32. A recent study33
city = 0.58, sensitivity = 0.67). No significant prediction showed that Ruminococcaceae was negatively correlated with the
power for the classification in the AED group (AUC = 0.53, fibrosis severity. In agreement with this, members in Rumino-
specificity = 0.40, sensitivity = 0.67) and the AEx group coccaceae was found negatively correlation with HFC in our
(AUC = 0.52, specificity = 0.50, sensitivity = 0.60) was found. study. In addition, some members of Bacteroides contribute to
Second, since ASVs abundance could not describe the individual the release of energy from dietary fiber and starch, and they are
network characteristics, we developed a Single SparCC network likely to be a major source of propionate34. Accordingly, we
method by combining a SparCC network26 with a sample specific found that members in Bacteroides was positively correlated with
network (SSN)30 (Fig. 5a). We used baseline metagenomics data propionic acid, isobutyric acid and isovaleric acid. A previous
to build a network for each patient (Fig. 5b) and found that study showed that compared to western countries, the abun-
responders in all three intervention groups tended to have more dance of Bacteroides in NAFLD is much lower in Chinese
interactions between species than low/non-responders, but no individuals35. However, our result suggested an important
statistical significance was observed (Fig. S5). Furthermore, in the energy-extracting role of members in Bacteroides, although their
linear regression analysis, the edge number in the AEx group abundance is relatively low in Chinese population. Of note, there
significantly predicted the change of HFC (Fig. 5c). Then we are currently conflicting results about the shifts of microbiota at
performed an unsupervised classifier to differentiate the respon- the taxon level in NAFLD-related studies. For example, some
ders from the low/non-responders by each Single SparCC net- studies showed that Bacteroides abundance was higher in NASH
work edge number. We found that AUC was about 70% in (nonalcoholic steatohepatitis, a stage of NAFLD progression)
different intervention groups (Fig. 5d), and the AUC of super- than in non-NASH36. Another article showed the phylum Bac-
vised Least absolute shrinkage and selection operator (LASSO) teroidetes (44.63%) tended to be more abundant in healthy
classifiers were higher in all three groups (Fig. 5e). For the pur- subjects than in NAFLD35. Our results showed that members in
pose of comparison, we also tested the ability of age, body weight Bacteroides abundance was increased after both diet and exercise
(WT) and body mass index (BMI) to distinguish the respon- intervention, and this increase was correlated with decreased
siveness of the interventions (Fig. S6). The result showed that HFC. However, previous studies showed higher fiber intake was
these clinical parameters were not able to differentiate the marginally associated with lower abundance of Bacteroides
responders from the low/non-responders (age: AUC = 0.437, uniformis37, and Bacteroides was increased in participants with
p = 0.347; WT: AUC = 0.470, p = 0.655; BMI: AUC = 0.519, obesity but decreased in lean participants after exercise
p = 0.771), except for BMI in the AED group. intervention38. These inconclusive results may be ascribed to
Considering the small sample size of metagenomic data, we differences in ethnicity, living environment and lifestyle among
performed the same Single SparCC network analysis on the the cohorts. Thus, the mechanism of action of the key microbiota
intervention groups’ samples by using the 16 S rRNA gene needs to be investigated in future studies.
sequencing data (Fig. S7). Consistent with the metagenomic The long-term adaptation after intervention may account for
results, the responders of the three intervention groups tended to maintenance of gut microbiota alpha diversity in all intervention
have more network edges than low/non-responders, but none of groups in contrast to the NI group in which alpha diversity
them were statistically significant (Fig. S8a). In regression decreased. However, the connectivity and robustness of
analysis, we found that, in addition to the AEx group, the the microbiota co-occurrence networks improved in three

NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications 7

ARTICLE NATURE COMMUNICATIONS | [Link]

a
1. Importing species abundance of N samples 2. Estimating Single SparCC network for each sample 3. Distinguishing responders, non/low-responders
Estimating the distribution of Single SparCC network by Single SparCC network
species abundance for each sample
1 2 N-1 N in N-1 samples Responders Non/low-responders

Species A 0 8 … 100 0
Species B 1 0 … 1000 0
Species C 100 2 … 8 30
Species A Species B
Estimating edges in Nth samples by comparing their
abundance to the N-1 distribution

b Responders Non-responders

AED

AEx

Diet

c d e
1.0

1.0
10 R2=0.26 P = 0.0914
R2=0.45 P = 0.0166
0.8

0 R2=0.35 P =0.0739 0.8

HFC change

0.6
Sensitivity

0.6
Sensitivity

AUC AUC
-20 AED 0.77 AED 0.90
AEx 0.74 AEx 0.94
0.4

Diet 0.72 0.4 Diet 0.76

-40
0.2

0.2

0.0
0.0

-60
3000 6000 9000 120000 1.0 0.8 0.6 0.4 0.2 0.0
edge number 1.0 0.8 0.6 0.4 0.2 0.0
Group AED AEx Diet Specificity
Specificity

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. b The 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. c Linear regressions of edge number
(Single SparCC network) with the change of HFC after intervention. AED: n = 12, AEx: n = 12, Diet: n = 10. d ROC 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. e ROC curve by
supervised Lasso model performance. AED: n = 12, AEx: n = 12, Diet: n = 10.

intervention groups. Various ecological microbiome studies39,40 lower functionality due to fewer taxa present that can support
have shown that healthy people have a higher average network their function in ecosystems39. In combination with these find-
degree and hence greater connectivity of their gut microbiota ings, our results indicate that microbiota as an ecosystem were
network41. A poorly developed microbial network usually has more healthy and stable after intervention. Interestingly, when

8 NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications

NATURE COMMUNICATIONS | [Link] ARTICLE

ASVs abundance is used to predict the change of HFC, it has the Comparing with the conventional machine learning methods,
best predictive power in the diet group, but is not predictive in the although the reference samples can be viewed as the training
exercise group. However, when assessed the robustness of co- samples, our method can further extract the second-order sta-
occurrence networks, the exercise intervention showed the tistical information (e.g., correlation or association or network) of
highest increment of robustness score. A possible explanation is each individual approximately, thus providing system-wide fea-
that diet intervention directly affects the structure and composi- tures of each test sample. In general, our results are important in
tion of gut microbiota through changes in nutrients, while exer- terms of how to assess the intervention effect on physiology and
cise intervention influences the interaction between bacteria pathophysiology and to develop personalized lifestyle treatment
through the regulation of host metabolism. In addition, this study for such patients through their individual gut microbiota
confirmed that co-occurrence networks of microbiota are scale- network.
free networks, and further defined four keystone taxa, which have Our study provides novel information on how exercise and diet
an prominently high degree in each network. A recent study intervention affect the gut microbiome for patients with comorbidity
found that keystone taxa that are highly connected in the conditions at both the population and individual levels. However,
microbiome can explain microbiome compositional turnover there are some limits in this study. Our sample size for metagenomics
better than all taxa combined42. Therefore, we speculate that what is relatively small but it is comparable with most of the studies44. The
exercise intervention improves not the gut microbiota diversity, dietary and exercise intervention methods are often different from
but the stability of the interaction network of the gut microbiota trial to trial, and even if the intervention methods are the same, the
ecosystem, thereby improving the metabolism of the host. The intervention contents and intensity are different which makes the
keystone taxa determined by the microbiota population network comparison with other trials challenging. In addition, since the trial
topology may be a new perspective beyond finding key bacteria was not originally designed to predict response efficiency with per-
by abundance and diversity comparison. sonalized microbiota network, our predictions were not validated by
No disease-specific drug/medical therapy for NAFLD is avail- independent samples. Nevertheless, this study provides a new per-
able due to its complexly pathogenetic mechanism, and man- spective in the research of gut microbiota, that is, in addition to
agement of lifestyle (such as exercise and diet) in patients is the focusing on the differential strains themselves, the topological prop-
only effcient invervention in clinical practices. However, studies erties of the microbial ecosystem community network are also of
have shown that not all participants respond to invervents of great significance. Using the network features of gut microbiota to
exercise or diet in a similar way16, which may be attributable to predict individual response to exercise and diet interventions may
the individual-variation and response of gut microbiota in help to choose more precise and effective treatment for different
patients. In some studies, 50% of the study participants show no NAFLD pateints.
change in their gut microbiota composition after an exercise
intervention17. Then prediction of inprovment of clinical phe- Methods
notypes in different patiens based on their gut microbiota would The study was approved by the Ethics Committee of Shanghai Institute of
be important to increase effcieny of NAFLD treatment. In the Nutrition (06.01.2013) and has been retrospectively registered in the International
Standard Randomized Controlled Trial Number Register (ISRCTN 42622771). We
previous studies, they ascribed to use the group-level changes of acknowledge that one limitation of the study is that the trial was restrospectively
bacteria abundance to construct the prediction models, which registered. Two of the primary outcomes (Glucose tolerance test and Liver fat
disregards the fact that the menbers of gut microbiota show content) have been previously reported21. All participants signed informed consent
ecological interactions such as competition or cooperation with forms before beginning the study. The study was conducted in accordance with the
each other. Our Single SparCC method evaluated not only principles of the Declaration of Helsinki.
microbiome composition and abundance, but also the interac-
tions between species. Such a network can be used to quantify the Subjects and design. This study was a randomized controlled trial. Detailed
information was given in our previous publication21. The enrollment of partici-
system-wide features, e.g., robustness or stability of the system, pants was between January 6 and December 25, 2013. The eligibility criteria for
and further predict intervention effects. Various studies28,43 have intervention groups were as follows: men or women aged 50–65 years with
adopted the concept of connectivity and robustness from network impaired fasting glucose (IFG between 5.6 to 6.9 mmol/L) or impaired glucose
theory to investigate the rule of microbial community in diseases tolerance (IGT between 7.8 to 11.0 mmol/L 2 hour after the intake of 75 g glucose),
and diagnosed as NAFLD by 1H MRS (liver fat > 5%) and by questionnaire that
from an ecological perspective. Actually, there is increasing evi- ongoing or recent alcohol consumption is <21 drinks on average per week in men
dence indicating a link between undesirable health conditions and <14 drinks on average per week in women; no chronic cardiovascular, serious
with altered microbial assembly process and a more fragile musculoskeletal or gastrointestinal problems and not on extreme diets; and for
microbiome network40,41. Our results indicated that the Single women, serum follicle-stimulating hormone level greater than 30 IU/L and last
menstruation more than 6 months ago but within 10 years. The exclusion criteria
SparCC network of gut microbiota is a valid method to predict included body mass index (BMI) > 38 kg/m2; serious cardiovascular or muscu-
the responders from low/non-responders for exercise interven- loskeletal problems; diagnosed with Type 1 diabetes and T2D; and mental illness.
tion and therefore opens a new avenue to assess the effect of In the initial study design, we estimated the sample size based on previous lit-
exercise intervention on gut microbiota at an individual level. erature of liver fat and microbiota as primary outcome variable21. We estimated
that 34 subjects in each group would have 85% power for mean comparison
Interestingly, we found that when differentially altered ASVs
between the groups and therefore we targeted to have 50 subjects in each group.
abundance is used to classify the responders and low/non- However, during the recruitment period, we found that only about 20% of the
responders, it has the best prediction power in the diet group but subjects met the inclusion criteria and due to the limited funding, we re-calculated
invalid in the exercise group. In terms of the AUC values, the the sample size based on our previous similar type of study45,46. Thus, for this
Single SparCC network of gut microbiota was better differential report, sample size estimation for the primary outcome HFC with 29 individuals
has 95% power to test against the null hypothesis that there is no change in any
features to responders from low/non-responders than age, body group. Further, when intervention groups compared with NI group having
weight and BMI in our study. These results indicate that corre- 17 subjects, the power for the HFC was 84% with < 0.05 two-sided significance
sponding gut microbial characteristics can be used to predict level. After every subject reached minimal 6 months intervention, the study did not
responsiveness of a subject to specific intervention methods. This continue to follow-up all the subjects and consequently, the follow-up study was
terminated. There were no sex distribution differences between the groups and the
is of great clinical significance in improving effectiveness and male/female ratio for each group was n = 6/23 for AEx, n = 6/22 for Diet, n = 7/22
efficiency of interventional therapy of NAFLD. Notably, the for AED, and n = 7/22 for NI group, respectively.
personalized prediction here requires a number of samples or
other individuals measured as the reference samples, against Interventions. The duration of the intervention ranged from 6.5 months to
which the network of each test sample can be projected. 11.1 months and on average 8.6 months. There were no significant differences in

NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications 9

ARTICLE NATURE COMMUNICATIONS | [Link]

the intervention duration among the groups (AEx 8.75 months, Diet 8.6 months, (Released 13.8) and then classified into ASVs at a threshold of 97% sequence
AED 8.5 months and NI 8.6 months) as well as between responders (8.6 months) identity using UCLUST. The representative sequence for each ASV was selected
and non-responders (8,7months). using default parameters and was imported into the latest Greengenes ARB
The exercise (AEx) group participated in a supervised progressive aerobic database to construct a phylogenetic tree.
exercise training program (such as Nordic brisk walking + stretching and other Raw sequencing data were analyzed by QIIME 2. In the process of running the
group exercises). Exercise was performed 2–3 times a week, 30–60 min per session DADA2 pipeline, based on the quality profile of the data, forward and reverse reads
with 60–75% of the maximum oxygen uptake (estimated from fitness test). Each were trimmed accordingly to ensure that the median quality score for each position
exercise session included a 5 min warm-up and 5 min cool-down period. The diet was above 32. The taxonomy of all ASVs were annotated by the SILVA (v132)
(Diet) group received a daily prepared meal (lunch), which accounted for 30–40% reference database47. All samples were rarefied to 10,000 per sample for
of the total daily energy intake on the basis of each individual’s dietary intakes and downstream analysis.
body weight. The meal included 37–40% carbohydrate with 9–13 g as fiber, 35–37%
fat (SAFA 10%, MUFA 15-20%, PUFA 19 10%) and 25–27% protein plus 5 g of Metagenomics analysis. Forty-two paired samples were sequenced at the
soluble fiber (dietary water-soluble fiber). The lunch box was then delivered to the
Novogene Bioinformatics Institute (Novogene, Beijing, China). After construction
study district community office where the study subjects were gathered. If the
of 350 bp insert libraries, metagenomic samples were pair-end sequenced on the
participants did not have time to come to pick up the meal, it was delivered to their
HiSeq X Ten System. Cluster generation, template hybridization, isothermal
home and the participants ate it for dinner. The other two meals (breakfast and
amplification, linearization, and blocking denaturing and hybridization of the
dinner) were cooked by subjects themselves following the nutritionist advice. In
sequencing primers were performed according to the workflow specified by the
addition, the Diet group was advised to maintain their current level of physical
service provider. The pipeline was consistent with previous studies48,49. We
activity during the intervention. The Exercise plus Diet (AED) group performed
removed contamination of host DNA and assembled the metagenome by using
the same exercise program and followed the same diet as described above for the
Bowtie (Default parameter Settings:—end-to-end,—sensitive, -i 200, -x 400) and
AEx and Diet groups. The no intervention (NI) group was advised to maintain
SOAP denovo (K-mers = 55)50.
their current level of physical activity and eating habits during the intervention.

Short chain fatty acid measurement. Slice of human feces were homogenized
Patient and public involvement. The patients were engaged via advertisements with 1.2 mL of phosphate buffer solution in weighed 2 mL centrifuge tubes (marker
distributed in 7 health clinical service centers in the Shanghai Yangpu district. In as W1). The mixture was vibrated completely and centrifuged twice at 4 °C at
addition, media and social media were used to inform about the study. Patients 16000 × g for 15 min. The supernatants were filtered and sterilized through a
were not involved in the design of the study. The study protocol was described to 0.22 μm nylon filter (EMD Millipore). Feces residue and tube were dried in 80 °C
study participants on the first study visit. drying apparatus overnight and weighed (marked as W2). 200 μL of the super-
natant was acidized with 100 μL of 50% (v/v) sulfuric acid. After vortexing and
Fecal DNA extraction. Faecal samples were collected at the baseline and after the standing for 2 min, the organic phase was extracted by adding 400 μL of diethyl
intervention. Microbial genome DNA was extracted from fecal samples using the ether. The extracted solution was measured by gas chromatography (GC) on an
bead-beating method (alterations in the gut microbiome and metabolism with Agilent 6890 (Agilent Technologies, CA, USA). The concentrations of 6 short chain
coronary artery disease severity). Briefly, we used an InviMag® Stool DNA kit fatty acids (acetic acid, propionic acid, butyric acid, isobutyric acid, valeric acid and
(Invitek, Berlin, Germany) with agitation in a mini-bead beater (Biospec Products, isovaleric acid) were calculated by each peak area and feces dry weight (W2-W1).
Bartlesville, OK, USA) as follows: approximately 0.2 g of thawed feces were added
to a 2 mL screw-cap tube containing 1 mL of lysis buffer P of the kit and 0.3 g Co-abundance gene groups (CAGs). The ASVs shared by at least 20% among all
Zirconia beads (0.1 mm, Biospec Products, Bartlesville, OK, USA). After sufficient samples were considered key ASVs. The correlation among 279 key ASVs was
homogenization by vortex for approximately 5 min, bead beating was performed calculated by the SparCC algorithm51 with a bootstrap procedure repeated 100
for 1 min at maximum speed. DNA was extracted by following the manufacturer’s times, and then correlation matrices were computed from the resampled data
instructions for bacterial DNA extraction involving proteinase K treatment and matrices. Once the bootstrapped correlation scores have been computed, only
subsequent purification using the KingFisher device (Invitek, Berlin, Germany). ASVs with correlation scores greater than 0.4 were classified into CAGs. The
p value at the desired cut-off < 0.05.
qPCR of total fecal bacteria. A plasmid of the 16 S full-length positive Rumi-
nococcus strain (n = 1010 copies/μl) was diluted according to different gradients Network fragility. Given the N node (CAG) number in a network, the N nodes
successively to 109, 108, 107, 106, 105, 104, 103, and 102 copies/μl. qPCR was were removed one by one in a random order. The fragility statistics (y axis) were
performed in a 20-μl reaction system containing template (20 ng), primer Uni331F the ratio of the node number in the largest connected component to the node
(5’-TCCTACGGGAGGCAGCAGT-3’), primer Uni797R(5′-GGACTACCAGGGT number N, and were calculated when each node was removed, while the percentage
A TCTAATCCTGTT-3′), and supermix (Bio-Rad) on a qTOWER3G touch sys- of removed nodes was taken as the x axis. The robustness value R for a network was
tem, with 2 replicates for standard and sample DNA. The PCR conditions were related to a function σ, which was dependent on the decreasing order of fragility
95 °C for 5 min, followed by 40 cycles of 95 °C for 20 s, 60 °C for 60 s, and 80 °C for statistics as nodes were removed28, i.e.,
5 s, and Melting curve 60 °C to 95 °C. A standard curve was determined though a
linear fifit of the copy number and CT value of the plasmid in different gradients.
1 N
The copy number of sample DNA was calculated through a standard curve. R¼ ∑ σði=NÞ
N i¼1
Construction of amplicon library of V3-V4 region of 16 S rRNA and sequen- where i represents the ith CAG. To ensure the results were stable, we repeatly
cing. The extracted DNA from each sample was used as the template to amplify the remove the nodes in each network for 1000 times. The median fragility statistics of
V3-V4 region of 16 s rRNA genes with bar-coded primers, forward primers-338F the 1000 curves were taken as the coordinates for the final curve to plot the
(5′- sampleIDtag- ACTCCTACGGGAGGCAGCA-3′) and reverse primers-806R robustness curve.
(5′- sampleIDtag-GGACTACHVGGGTWTCTAAT). PCR reactions were run in a
thermocycler PCR system (Applied Biosystems, USA) using the following pro-
gramme: 5 min of denaturation at 96 °C followed by 25 cycles of for 30 sec at 96 °C Single SparCC network analysis. The traditional correlation coefficient mainly
(denaturation), 30 sec at 50 °C (annealing) and 30 sec at 72 °C (elongation), with a exploits the general information in a group of samples. Rather than focus on a
final extension at 72 °C for 5 min. PCR product was excised from a 1.5% agarose gel ‘mean’ value within a group, our Single SparCC network is proposed to investigate
and purified by QIAquick Gel Extraction Kit (QIAGEN, cat# 28706). Purified PCR the characteristics of each sample. Gut microbiota 16 S sequencing and metage-
products from the fecal samples were combined at equal concentrations and used nomics data suffered from spurious correlation. SparCC is the most popular
to construct a metagenomic library using an Illumina TruSeq sample preparation method to eliminate spurious correlation26. SparCC can detect the association
kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. between species in a set of sample size > 1, but it is invalid for an individual sample.
Barcoded V3-V4 amplicons were sequenced using a pair-end method via the Our approach enabled us to explore the personalized characteristics of NAFLD
Illumina Miseq (Illumina, San Diego, CA, USA) sequencing platform with a 6-cycle patients.
index read. The correlation of species i and j can be inferred from their ratio. Even though
high-throughput sequencing only produce relative quantification of microbe, the
ratio of those relative value is the same with real abundance. Define tij to describe
16 S rRNA Sequence analysis. Quality control of the raw data was performed as variance of log transformed relative value of component xi and xj:
follows: (1) the sequence should have a perfect match to the barcode in at least one
end; (2) the sequence should have a BLAST match to at least one end of the 16 S " # " #
ωi xi
rRNA gene V3-V4 region primers; (3) the length of the trimmed sequence (without t ij  Var log ¼ Var log ; ð1Þ
barcodes/primers) should be between 400 nt and 500 nt; and (4) there should be no ωj xj
more than two undetermined bases. All the retained sequences were processed
using the QIIME 2 (Quantitative Insights into Microbial Ecology) package. The where xi and xj are the relative abundances of microbe x and y, while ωi and ωj are
sequences were aligned using the PyNAST aligner with the Greengenes core set the absolute abundances of those two microbes.

10 NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications

NATURE COMMUNICATIONS | [Link] ARTICLE

In SparCC method, tij was related to Pearson correlation coefficient ρij by community membership (weighted UniFrac distance for bacterial relative abun-
equation: dance) and functional composition (bray curtis distances for KEGG pathway).
Pairwise comparisons of beta diversity between groups were performed using
t ij ¼ ω2i þ ω2j  2ρij ωi ωj ; ð2Þ [Link] (R package RVAideMemoire v0.9). p < 0.05 is considered as
significant in all the above analyses. The Spearman correlation analysis between
where ω2i and ω2j are the variances of the log-transformed real abundance of species i
microbiome and clinical parameters was performed with MATLAB R2020a in
and j. Because of the sparse nature of microbial association network, we can assume samples with complete data. Co-occurrence network analysis was conducted for
that the average of all ρij correlation is small enough to ignore, which allow us to solve samples before and after intervention in all 4 groups using pairwise Spearman
for ω2i and ω2j . Then put in the values of ω2i and ω2j back to Eq. (2), and end up with correlations in MATLAB R2020a, and only the significant correlations (p < 0.05)
the approximated correlation ρij. Suppose the sample size is S and the number of nodes were used for network construction. The networks in each intervention group
is D. Letting a = 2D-3, b ¼ 2ðD1ÞðD2Þ 1
, we achieve an analytical solution for ρij: (before and after intervention) were merged together to facilitate direct comparison
!2 and were further visualized in Cytoscape software (v3.8.1). To detect differential
1 S D xðkÞ taxonomy, KEGG pathway analysis and linear discriminant analysis effect size
ρij ¼ ∑ bða þ 1Þ ∑ log iðkÞ  μlogxij
2sωi ωj k¼1 j¼1 xj (LEfSe v1.0, which is an algorithm for high-dimensional biomarker discovery and
!2 !2 identifies genomic features such as genes, pathways or taxa as well as characterizing
ðkÞ
D D x D xðkÞ the differences between two or more biological conditions or classes) was per-
 2b ∑ ∑ log iðkÞ  μlogxij þ bða þ 1Þ ∑ log iðkÞ  μlogxij ð3Þ formed with cut-offs (LDA Score > 2.0, p < 0.05).
j¼1 i¼1 xj i¼1 xj
!2 !!
xðkÞ Reporting summary. Further information on research design is available in the Nature
 log iðkÞ  μlogxij
xj Research Reporting Summary linked to this article.

where
!2 Data availability
1 D S xði kÞ The study protocol has been published previously ([Link]
ω2i ¼ bða þ 1Þ  ∑ ∑ log ðkÞ  μlogxij 14-48). The raw Illumina sequence data including 16 s rRNA sequence and
S j¼1 k¼1 xj
ð4Þ metagenomics generated in this study have been deposited to the NCBI database under
!2
1 D D S x kÞ
ð BioProject accession no. PRJNA757939 ([Link]
 b  ∑ ∑ ∑ log iðkÞ  μlogxij : term=PRJNA757939). The data is available for academic use under controlled access in
S i¼1 j¼1 k¼1 xj compliance with the regulation of the Ministry of Science and Technology (MOST) of
China for the deposit and use of human genomic data. All individual de-identified
!
D S xði kÞ participant data that support the findings of this study are available from the
μlogxij ¼ ∑ ∑ log ; ð5Þ corresponding author upon reasonable request based on the rule of data protection and
j¼1 k¼1 xðj kÞ consent. The access to the controlled data will be valid for one year from the time of the
Since ρij is a form of k elements summation, we can break it down into a vector data accessibility approved. The processed data are available within the Source Data file.
that depends on the correlation coefficient ρij, and each element ρðkÞ
ij is the sample
specific correlation: Code availability
!2 !2 ! Codes used for computing population and individual networks are available on GitHub
1 D xði kÞ D D xði kÞ
ρðkÞ ¼ bða þ 1Þ ∑ log  μlogx  2b ∑ ∑ log  μlogx ([Link]
2sωi ωj xðj kÞ xðj kÞ
ij
j¼1 ij j¼1 i¼1 ij
!2 !2 ! ð6Þ
D xði kÞ xði kÞ
þ bða þ 1Þ ∑ log  μlogx  log  μlogx : Received: 21 February 2021; Accepted: 11 April 2022;
i¼1 xðj kÞ ij xðj kÞ ij

Thus, the algorithm is given as follows:

Step 1: From the taxonomy abundance of s samples, we calculated μlogxij , ωi and ωj.
Step 2: To infer the correlation ρðkÞ
ij in the network for kth sample, we substitute
μlogxij , ωi and ωj into Eq. (6).
Step 3: We iteratively calculate the correlation coefficients for each sample. References
1. Rinella, M. E. Nonalcoholic fatty liver disease: a systematic review. JAMA 313,
2263–2273 (2015).
Prediction of responders in change of HFC. In the supervised prediction process,
2. Williams, C. D. et al. Prevalence of nonalcoholic fatty liver disease and
we used baseline samples from each group to predict the responders from low/non-
responders after intervention. The analysis takes three steps. nonalcoholic steatohepatitis among a largely middle-aged population utilizing
Step 1: Through the Single SparCC algorithm, we generated one SSN for each ultrasound and liver biopsy: a prospective study. Gastroenterology 140,
individual in the n samples. We then converted these networks into an edge matrix, 124–131 (2011).
which contained the variables consisting of 3321 edges. The edge weight in n 3. Leite, N. C., Salles, G. F., Araujo, A. L., Villela-Nogueira, C. A. & Cardoso, C.
samples was taken as the feature. R. Prevalence and associated factors of non-alcoholic fatty liver disease in
Step 2: We used the N-sample edge matrix as the input set. The level of HFC patients with type-2 diabetes mellitus. Liver Int. 29, 113–119 (2009).
after the intervention was the dependent variable, which was used to train an 4. Lonardo, A., Ballestri, S., Marchesini, G., Angulo, P. & Loria, P. Nonalcoholic
Elastic-Net model. fatty liver disease: a precursor of the metabolic syndrome. Dig. Liver Dis. 47,
Step 3: Based on the predicted and true values from step 2, we calculated the 181–190 (2015).
prediction metrics. Thereafter, responders from low/non-responders according to 5. Betrapally, N. S., Gillevet, P. M. & Bajaj, J. S. Gut microbiome and liver
their change of HFC after intervention were identified. disease. Transl. Res 179, 49–59 (2017).
In addition, we performed a regression analysis to use the Single SparCC 6. Qin, N. et al. Alterations of the human gut microbiome in liver cirrhosis.
network edges to assess whether they predict HFC change after intervention by Nature 513, 59–64 (2014).
using the stat_poly_eq function in the R package ‘ggpmisc’. 7. Betrapally, N. S., Gillevet, P. M. & Bajaj, J. S. Changes in the intestinal
In an unsupervised prediction process, we directly used the network attributes microbiome and alcoholic and nonalcoholic liver diseases: causes or effects?
(edge number, mean degree) of a Single SparCC network (before intervention) to Gastroenterology 150, 1745–1755.e1743 (2016).
predict the responders from the low/non-responders. A ROC curve was then 8. Haas, J. T., Francque, S. & Staels, B. Pathophysiology and mechanisms of
drawn to visualize the predictive effect (roc function in R package pROC). nonalcoholic fatty liver disease. Annu. Rev. Physiol. 78, 181–205 (2016).
9. Ratziu, V., Bellentani, S., Cortez-Pinto, H., Day, C. & Marchesini, G. A
Statistical analysis. Before statistical test for relative abundance data, taxonomy position statement on NAFLD/NASH based on the EASL 2009 special
that appeared in less than 20% of the samples were excluded. Alpha diversity conference. J. Hepatol. 53, 372–384 (2010).
(Shannon index) was performed using the phyloseq (v1.33) package in R v4.0. For 10. Mardinoglu, A. et al. An integrated understanding of the rapid metabolic
the comparison of Shannon index and the absolute content of microbiome, general benefits of a carbohydrate-restricted diet on hepatic steatosis in humans. Cell
linear model of analysis of variance controlled for covariates (ANCOVA) for Metab. 27, 559–571.e555 (2018).
repeated measures (2 factor interactions: group x time) and controlled for change 11. Johnson, N. A. et al. Aerobic exercise training reduces hepatic and visceral lipids
of body weight, baseline value and intervention duration followed by Sidak cor- in obese individuals without weight loss. Hepatology 50, 1105–1112 (2009).
rection for multiple comparison between the groups. Contrast results (K Matrix) 12. Oh, S. et al. Moderate to vigorous physical activity volume is an important
were used to localize the group differences: *p < 0.05, **p < 0.01 and ***p < 0.001 factor for managing nonalcoholic fatty liver disease: a retrospective study.
by contract to the NI group. Beta diversity analyses were performed to assess the Hepatology 61, 1205–1215 (2015).

NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications 11

ARTICLE NATURE COMMUNICATIONS | [Link]

13. Clarke, S. F. et al. Exercise and associated dietary extremes impact on gut insomnia symptoms: a randomized controlled trial. Sleep. Med 25, 113–121
microbial diversity. Gut 63, 1913–1920 (2014). (2016).
14. Mailing, L. J., Allen, J. M., Buford, T. W., Fields, C. J. & Woods, J. A. Exercise 46. Tan, X. et al. Effect of six-month diet intervention on sleep among overweight
and the gut microbiome: a review of the evidence, potential mechanisms, and and obese men with chronic insomnia symptoms: a randomized controlled
implications for human health. Exerc Sport Sci. Rev. 47, 75–85 (2019). trial. Nutrients 8, (2016).
15. Codella, R., Luzi, L. & Terruzzi, I. Exercise has the guts: how physical activity 47. Quast, C. et al. The SILVA ribosomal RNA gene database project: improved
may positively modulate gut microbiota in chronic and immune-based data processing and web-based tools. Nucleic Acids Res. 41, D590–596 (2013).
diseases. Dig. Liver Dis. 50, 331–341 (2018). 48. Xing, J. et al. Hypoxia induces senescence of bone marrow mesenchymal stem
16. Anjo, D. et al. Who are the non-responder patients to cardiac rehabilitation? cells via altered gut microbiota. Nat. Commun. 9, 2020–2020 (2018).
Eur. Heart J. 34, [Link] (2013). 49. Li, J. et al. Gut microbiota dysbiosis contributes to the development of
17. Mohr, A. E. et al. The athletic gut microbiota. J. Int. Soc. Sports Nutr. 17, 24 hypertension. Microbiome 5, 14 (2017).
(2020). 50. Luo, R. et al. SOAPdenovo2: an empirically improved memory-efficient short-
18. Kirwan, L. et al. Phenotypic factors influencing the variation in response of read de novo assembler. GigaScience 1, 18 (2012).
circulating cholesterol level to personalised dietary advice in the Food4Me 51. Wishart, D. S. et al. HMDB 4.0: the human metabolome database for 2018.
study. Br. J. Nutr. 116, 2011–2019 (2016). Nucleic Acids Res. 46, D608–d617 (2018).
19. Sparks, L. M. Exercise training response heterogeneity: physiological and
molecular insights. Diabetologia 60, 2329–2336 (2017).
20. Banerjee, S., Schlaeppi, K. & van der Heijden, M. G. A. Keystone taxa as Acknowledgements
drivers of microbiome structure and functioning. Nat. Rev. Microbiol. 16, Funding to S.C. was provided by the China State Sport General Administration (2013B040,
567–576 (2018). 2015B039), the National Natural Science Foundation of China (Nos. 31571219), the
21. Cheng, S. et al. Effect of aerobic exercise and diet on liver fat in pre-diabetic Shanghai Jiao Tong University Zhiyuan Foundation (CP2014013), and to S.L. the Inter-
patients with non-alcoholic-fatty-liver-disease: A randomized controlled trial. disciplinary Program of Shanghai Jiao Tong University (ZH2018QNB05). Funding to
Sci. Rep. 7, 15952 (2017). L.C. was supported by the National Key R&D Program of China (No. 2017YFA0505500),
22. Segata, N. et al. Metagenomic biomarker discovery and explanation. Genome Strategic Priority Research Program of the Chinese Academy of Sciences (No.
Biol. 12, R60 (2011). XDB38040400), National Natural Science Foundation of China (Nos. 31930022, 31771476,
23. Koh, A., De Vadder, F., Kovatcheva-Datchary, P. & Bäckhed, F. From dietary 12026608). Funding to C.Z. was supported by the National Natural Science Foundation of
fiber to host physiology: short-chain fatty acids as key bacterial metabolites. China (Nos. 31930022 and 81871091).In addition to the authors, we thank Mengya Lin for
Cell 165, 1332–1345 (2016). contribution to the absolute microbial abundance measurement. We would like to thank the
24. Nielsen, H. B. et al. Identification and assembly of genomes and genetic study team members and participants for their contributions to the success of this trial.
elements in complex metagenomic samples without using reference genomes.
Nat. Biotechnol. 32, 822–828 (2014). Author contributions
25. Zhang, C. et al. Dietary modulation of gut microbiota contributes to S.C. designed the study and has full access to all of the data in the study and takes full
alleviation of both genetic and simple obesity in children. EBioMedicine 2, responsibility for the integrity of the data and for the accuracy of the data analysis. S.C., C.Z.,
968–984 (2015). S.L., J.G., D.L. participated in subject recruitment, measurements and intervention process.
26. Friedman, J. & Alm, E. J. Inferring correlation networks from genomic survey Y.Y., R.C., S.L. and X.Y. processed fecal samples and the 16 s rRNA sequence analysis. L.W.
data. PLoS Comput Biol. 8, e1002687 (2012). and L.C. involved in the metagenomics data analysis. X.Z. performed SCFA analysis. R.C.,
27. Faust, K. & Raes, J. Microbial interactions: from networks to models. Nat. Rev. L.W., Y.Y. and T.X. performed the statistical analyses. L.X., P.W., C.Z., L.C. and S.C. checked
Microbiol 10, 538–550 (2012). data analyses and interpretation of the results. R.C., L.W. and S.C. drafted the manuscript.
28. Iyer, S., Killingback, T., Sundaram, B. & Wang, Z. Attack robustness and L.X., P.W., C.Z., L.C. and S.C. performed critical revision of the manuscript. All authors have
centrality of complex networks. PLoS ONE 8, e59613 (2013). approved the final submission version of the manuscript.
29. Strogatz, S. H. Exploring complex networks. Nature 410, 268–276 (2001).
30. Liu, X., Wang, Y., Ji, H., Aihara, K. & Chen, L. Personalized characterization
of diseases using sample-specific networks. Nucleic Acids Res. 44, e164 (2016). Competing interests
31. Christopherson, M. R. et al. Unique aspects of fiber degradation by the The authors declared no conflict of interest.
ruminal ethanologen Ruminococcus albus 7 revealed by physiological and
transcriptomic analysis. BMC Genomics 15, 1066 (2014).
32. Crost, E. H. et al. Utilisation of mucin glycans by the human gut symbiont
Additional information
Supplementary information The online version contains supplementary material
Ruminococcus gnavus is strain-dependent. PLoS ONE 8, e76341 (2013).
available at [Link]
33. Lee, G. et al. Distinct signatures of gut microbiome and metabolites associated
with significant fibrosis in non-obese NAFLD. Nat. Commun. 11, 4982 (2020).
Correspondence and requests for materials should be addressed to Chenhong Zhang,
34. Flint, H. J. & Duncan, S. H. In Encyclopedia of Food Microbiology (Second
Luonan Chen or Sulin Cheng.
Edition) (eds Carl A. Batt & Mary Lou Tortorello) 203–208 (Academic Press,
2014).
Peer review information Nature Communications thanks Yinghui Wei and the other
35. Shen, F. et al. Gut microbiota dysbiosis in patients with non-alcoholic fatty
anonymous reviewer(s) for their contribution to the peer review of this work. Peer
liver disease. Hepatobiliary Pancreat. Dis. Int. 16, 375–381 (2017).
reviewer reports are available.
36. Boursier, J. et al. The severity of nonalcoholic fatty liver disease is associated
with gut dysbiosis and shift in the metabolic function of the gut microbiota.
Reprints and permission information is available at [Link]
Hepatology 63, 764–775 (2016).
37. Lin, D. et al. Association of dietary fibre intake and gut microbiota in adults. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in
Br. J. Nutr. 120, 1014–1022 (2018). published maps and institutional affiliations.
38. Allen, J. M. et al. Exercise alters gut microbiota composition and function in
lean and obese humans. Med Sci. Sports Exerc 50, 747–757 (2018).
39. Wagg, C., Schlaeppi, K., Banerjee, S., Kuramae, E. E. & van der Heijden, M. G.
A. Fungal-bacterial diversity and microbiome complexity predict ecosystem Open Access This article is licensed under a Creative Commons
functioning. Nat. Commun. 10, 4841 (2019). Attribution 4.0 International License, which permits use, sharing,
40. Strona, G. & Lafferty, K. D. Environmental change makes robust ecological adaptation, distribution and reproduction in any medium or format, as long as you give
networks fragile. Nat. Commun. 7, 12462 (2016). appropriate credit to the original author(s) and the source, provide a link to the Creative
41. Kim, H.-J. et al. Fragile skin microbiomes in megacities are assembled by a Commons license, and indicate if changes were made. The images or other third party
predominantly niche-based process. Sci. Adv. 4, e1701581 (2018). material in this article are included in the article’s Creative Commons license, unless
42. Herren, C. M. & McMahon, K. D. Keystone taxa predict compositional change indicated otherwise in a credit line to the material. If material is not included in the
in microbial communities. Environ. Microbiol 20, 2207–2217 (2018). article’s Creative Commons license and your intended use is not permitted by statutory
43. Kurtz, Z. D. et al. Sparse and compositionally robust inference of microbial regulation or exceeds the permitted use, you will need to obtain permission directly from
ecological networks. PLoS Comput Biol. 11, e1004226 (2015). the copyright holder. To view a copy of this license, visit [Link]
44. Liu, Y. et al. Gut microbiome fermentation determines the efficacy of exercise licenses/by/4.0/.
for diabetes prevention. Cell Metab. 31, 77–91.e75 (2020).
45. Tan, X., Alén, M., Wiklund, P., Partinen, M. & Cheng, S. Effects of aerobic
exercise on home-based sleep among overweight and obese men with chronic © The Author(s) 2022

12 NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications

NATURE COMMUNICATIONS | [Link] ARTICLE

1
Exercise Translational Medicine Center, Shanghai Center for Systems Biomedicine, Shanghai Jiao Tong University, Shanghai, China. 2Exercise,
Health and Technology Center, Faculty of Physical Education, Shanghai Jiao Tong University, Shanghai, China. 3School of Life Sciences and
Biotechnology, Shanghai Jiao Tong University, Shanghai, China. 4Key Laboratory of Systems Biology, Shanghai Institute of Biochemistry and Cell
Biology, Center for Excellence in Molecular Cell Science, Chinese Academy of Sciences, Shanghai, China. 5Faculty of Sport Sciences, University of
Jyväskylä, Jyväskylä, Finland. 6School of Physical Education and Training, Shanghai University of Sport, Shanghai, China. 7Ningbo University,
School of Medicine, Ningbo, China. 8Shidong Hospital of Yangpu District, Shanghai, China. 9School of Kinesiology, Shanghai University of Sport,
Shanghai, China. 10State Key Laboratory of Microbial Metabolism, School of Life Sciences and Biotechnology, Shanghai Jiao Tong University,
Shanghai, China. 11Key Laboratory of Systems Health Science of Zhejiang Province, Hangzhou Institute for Advanced Study, University of Chinese
Academy of Sciences, Hangzhou, China. 12Guangdong Institute of Intelligence Science and Technology, Zhuhai, China. 13These authors
contributed equally: Runtan Cheng, Lu Wang. ✉email: zhangchenhong@[Link]; lnchen@[Link]; [Link]@jyu.fi

NATURE COMMUNICATIONS | (2022)13:2555 | [Link] | [Link]/naturecommunications 13