Changes of gut microbiome composition and metabolites associated with hypertensive heart failure rats

The potential role of the gut microbiome (GM) in heart failure (HF) had recently been revealed. However, the underlying mechanisms of the GM and fecal metabolome in HF have not been characterized. The Dahl salt-sensitive rat model of hypertensive heart failure (H-HF) was used to study the clinical symptoms and characteristics. To elucidate the pathogenesis of HF, we combined 16S rRNA gene sequencing and metabolomics to analyze gut microbial compositions and fecal metabolomic profiles of rats with H-HF. PCoA of beta diversity shown that the gut microbiome composition profiles among the three groups were separated. Gut microbial composition was significantly altered in H-HF rats, the ratio of Firmicutes to Bacteroidetes(F/B) increased and the abundance of Muribaculaceae, Lachnospiraceae, and Lactobacillaceae decreased. Significantly altered levels of 17 genera and 35 metabolites were identified as the potential biomarker of H-HF. Correlation analysis revealed that specific altered genera were strongly correlated with changed fecal metabolites. The reduction in short-chain fatty acids (SCFA)-producing bacteria and trimethylamine N-oxide (TMAO) might be a notable characteristic for H-HF. This is the first study to characterize the fecal microbiome of hypertensive heart failure by integrating 16S rRNA gene sequencing and LC–MS-based metabolomics approaches. Collectively, the results suggesting changes of gut microbiome composition and metabolites are associated with hypertensive heart failure rats.


Background
Heart failure (HF) is the terminal stage of all cardiac diseases, with high morbidity and mortality rates [1]. Epidemiological data revealed that the prevalence of HF is 1-2% in adults and increases to more than 10% in people over the age of 70 [2]. The leading causes of HF are hypertension and ischemic heart diseases, and HF resulting from hypertension has become to be a major public health concern [3]. Therefore, timely diagnosis and early treatment are the keys. The gut microbiome (GM) is a complex community of trillions of bacteria in the gastrointestinal tract and has emerged as a central factor affecting human health and disease [4]. Hostmicrobiota interactions involving inflammatory and metabolic pathways have been linked to the pathogenesis of cardiovascular disease [5,6]. A growing number of studies have shown that GM is closely related to the occurrence and development of HF, in addition to alterations in GM composition, the metabolic potential of GM has been identified as a contributing factor in the development of diseases, including the trimethylamine (TMA)/ trimethylamine N-oxide (TMAO) pathway, short-chain fatty acids (SCFA) pathway, bile acid pathway and uremic toxin pathways, so microbiota is expected to become an essential target for intervention of HF [7]. HF has long been recognized to be associated with changed in intestinal function [8]. The gut hypothesis suggests that HF can lead to bacterial translocation across the intestine, and increase the level of bacteria throughout the systemic circulation with increased inflammatory status, thus promoting the further development of HF [9]. These studies had significantly increased attention towards the connection between our gut and heart. Thus, recognition of the gut-heart axis may lead to new insights and breakthroughs in the therapies for HF [10,11].
However, comprehensive analysis of the composition and metabolism of the microbiome in HF has not been conducted. Thus, studies are needed to investigate the fecal microbiome in association with HF and further reveal the effects of fecal metabolic changes in disease pathogenesis. Dahl salt-sensitive rat model is a wellestablished model of hypertensive heart failure [12][13][14]. Therefore, in the present study, we performed animal studies using Dahl salt-sensitive rats to evaluate intestinal microbial communities and metabolic profiles of H-HF, using 16S rRNA gene sequencing and metabolomics, to clarify the pathogenesis and consequences of HF.

Result
Echocardiographic and blood pressure measurement The echocardiographic parameters were shown in Fig. 1a and b,compared to the CON group and the SR group, both LVEF and LVFS decreased in H-HF group, suggesting compromised cardiac function, consistent with extant literature [14]. SBP and DBP were significantly higher in H-HF group, whereas blood pressure in the Fig. 1 Echocardiographic data on LVEF and LVFS, serum concentration of Nt-proBNP, Zonulin, LPS, and blood pressure. a: left ventricular ejection fraction (LVEF); b: left ventricular fractional shortening (LVFS). c: systolic blood pressure (SBP); d: diastolic blood pressure (DBP); e: NT-proBNP was assessed by ELISA; f: LPS was assessed by ELISA. g: Zonulin was assessed by ELISA. # P < 0.05, ## P < 0.01 compared with the CON group; *P < 0.05,** P < 0.01 compared with the SR group; Data are presented as the mean ± SD; n = 8 CON group and the SR group did not significantly change over time ( Fig. 1c and d).

Concentration of NT-proBNP, LPS, Zonulin in serum
The N-terminal proB-type natriuretic peptide (NT-proBNP) concentration is a sensitive and reliable biomarker for the diagnosis of heart failure. Our results revealed a significant increase in the NT-proBNP in the H-HF group compared to the CON group and the SR group (Fig. 1e). The intestinal barrier function is crucial for gut homeostasis. Serum lipopolysaccharide (LPS) was measured as an indicator of intestinal barrier function, and Zonulin was considered a marker of intestinal permeability. As shown in Fig. 1f and g, compared to the CON group and the SR group, the LPS and the Zonulin in the H-HF group significantly increased, indicating increased intestinal permeability and compromised intestinal barrier in the H-HF group rats.

Histological results
HE staining revealed the distribution of myocardial cells was regular in both CON group and SR group, while the myocardial cells in the H-HF group were swollen with irregular shapes and disordered arrangements, and substantial inflammatory cellular infiltrate in the myocardial interstitium was observed (Fig. 2a-c). HE-stained colonic section of the H-HF group exhibited integrity loss of the intestinal mucosa and infiltration of inflammatory cells into the colon tissue with lymphoid hyperplasia, which was not observed in the CON and SR groups ( Fig. 2 d-f). This suggests intestinal barrier damage (leaky gut) in the H-HF group.

Changes in metabolomics features of H-HF
PCA was applied to visualize metabolic alterations of the three experimental groups. As shown in the Fig. 3a, the CON group, H-HF group and SR group can be clearly separated. A segregation was visible between the CON and H-HF groups, indicating that certain significant biochemical changes occurred after the high salt diet. In addition, further OPLS-DA analysis also displays that the three groups were obviously separated (Fig. 3c,e). Meanwhile, a statistical validation of the OPLS-DA model was performed using 200 permutation tests, as shown in Fig. 3d and f, the model was reliable.
The metabolites with VIP > 1.0 and p < 0.05 were considered as significantly changed. In comparison of the H-HF group and CON group, 70 metabolites showed significantly different levels (Table S1), and 97 metabolites in the SR group were significantly different from those in the H-HF group (Table S2). The common differential metabolites were analyzed by Venn diagram (Fig. 3b), and it revealed that 35 metabolites were common of the three groups, the results are shown in Table 1. These metabolites could be regarded as potential biomarkers of H-HF. The heatmap of differential metabolite biomarkers is displayed in Fig. 4e. Differentially expressed biomarkers mainly consisted of organonitrogen compounds, organoheterocyclic compounds, organic oxygen compounds, organic nitrogen compounds, organic acids and derivatives, lipids and lipidlike molecules, and alkaloids and derivatives. Among them, studies have shown that metabolites associated with gut microbiota metabolism included creatinine, lithocholic acid, cholic acid, capric acid, glutaric acid, trimethylamine N-oxide, choline, and betaine [15]. Key metabolic pathway analysis for different metabolites All 35 potential biomarkers were subjected to metabolic pathway analysis (MetPA) using the KEGG online database and MetaboAnalyst. As shown in Fig. 4a, five metabolic pathways were influenced (P < 0.05), involving: (1) Histidine metabolism, (2) Arginine and proline metabolism, (3) Alanine, aspartate and glutamate metabolism, (4) Glycine, serine and threonine metabolism, and (5) Glycerophospholipid metabolism. The information for all the affected pathways is described in Table S3.

Diversity analysis of gut microbiota
A total of 1,920,409 paired-end reads were obtained from the 24 samples after sequencing (Table S4), with 1, 795,548 clean tags after alignment and filtering. Samples contained 74,814 clean tags on average. A total of 609 OTUs (gamma diversity) were obtained at a similarity level of 97%. The community composition of the three groups at the phylum (Fig. 5a) and family (Fig. 5b) was determined. For simplicity, only the 10 most abundant taxa are shown, with all others grouped into "other." Compared with the CON group and SR group, the abundance of Firmicutes in H-HF group increased significantly ( Fig. 5c), while Bacteroidetes decreased. Thus, the ratio of Firmicutes to Bacteroidetes increased significantly (Fig. 5d). An increase in the abundance of Ruminococcaceae, contrary to a decrease in the abundance of Muribaculaceae, Lachnospiraceae, and Lactobacillaceae, was observed in the H-HF group compared with the CON group and SR group (Fig. 5e). It revealed that the structure of gut microbiota in H-HF model has changed. Research suggests that the imbalance in the F/B ratio is largely related to energy metabolism [16,17], an increase in the F/B ratio in our study also indicates the Fig. 3 Results of multivariate analysis among three groups. Each dot with three kinds of color (green, H-HF model group; red, CON group; purple, SR group) represented the different sample. a: PCA score plots of among three groups (R 2 X = 0.632); b: Venn diagram showed the overlapping and unique differential metabolites amongst the comparison groups. c: OPLS-DA score plots of CON and H-HF group (R 2 X = 0.505, R 2 Y = 0.997, Q 2 = 0.977); d:The permutation test(n = 200) for the OPLS-DA model of CON and H-HF group; e: OPLS-DA score plots of H-HF and SR group (R 2 X = 0.387, R 2 Y = 0.976, Q 2 = 0.736); f: The permutation test (n = 200) for the OPLS-DA model of H-HF and SR group dysregulated energy metabolism, indicating that abnormal energy metabolism was found in H-HF [18], consistent with the existing literature.
Rarefaction curves (Fig. 5f) revealed the amount of sequencing data was sufficient to reflect the true taxonomic diversity. The value of Good's coverage for each group was higher than 99.8% (Table S4). To investigate the variances of gut microbiota structural diversity in the three groups, Chao1 and Shannon indexes of alpha diversity analysis were used to evaluate richness, evenness and diversity, respectively. Alpha diversity analysis revealed no significant difference in gut microbiota diversity between groups based on Chao 1 and Shannon indices (Fig. S1b,c), but unweighted unifrac PCoA of beta diversity shown that the gut microbiome composition profiles among the SR, CON and H-HF groups were separated (Fig. 4b). ANOSIM tests on UniFrac distance data were applied to permutational multivariate analysis(R = 0.753, P < 0.01). That is, the structural diversity of the gut microbiota was significantly different in H-HF group. Our results indicated a significant gut microbial shift during the development of HF, in accord with literatures [19], which highlights the potential role of gut microbiota in HF pathogenesis.

The alterations of the gut microbiota in H-HF
A total of 61 genera were found to vary significantly different in the H-HF and CON groups (Table S5) and 27 genera were found to differ be significantly different between in the SR and H-HF groups (Table S6). Venn diagram showed the overlapping and unique differential flora amongst the comparison groups. As shown in the Fig. 4c, there were a total of 17 different flora among the three groups. The microflora can be considered as the Fig. 4 a: Summary of the pathway analysis with MetPA of differential metabolites; b: PCoA analysis of three groups; c: Venn diagrams of three groups; d: Heat map of the differential gut microbiota of three groups; e: Heat map of the differential metabolites of three groups biomarkers of hypertensive heart failure. Variations in the identified differential gut microbiota of three groups were depicted in the heatmap (Fig. 4d).

Correlation of the gut microbiota and fecal metabolic
Pearson correlation analysis was performed for the screened fecal metabolites and the gut microbiome taxa at genus levels. The significantly related metabolites and gut genera are shown in the form of heat maps (Fig. 6). The related dataset of Pearson's correlation is described in Table S7. The Pearson's r-value > 0.75 is considered as a strong correlation. As shown in Fig. 7, for example, Fig. 5 The histogram of species distribution at the phylum (a) and family (b) levels in three groups; c:The relative abundance of Firmicutes and Bacteroidetes; d The F/B ratio; e: The relative abundance of Ruminococcaceae, Muribaculaceae, Lachnospiraceae, and Lactobacillaceae, f: Rarefaction curves;; # P < 0.05, ## P < 0.01 compared with the CON group;*P < 0.05, ** P < 0.01 compared with the SR group; data are presented as the mean ± SD; n = 8 Pyruvic acid and D-Maltose were strong negatively related to Adlercreutzia, Gordonibacter and uncultured_ bacterium_f_Erysipelotrichaceae, while Creatinine displayed strong negative correlation with Gordonibacter, and uncultured_bacterium_f_Erysipelotrichaceae. Pyruvic acid, D-Maltose and Creatinine are important indicator of energy metabolism. It indicates that Adlercreutzia, Gordonibacter and uncultured_bacterium_f_Erysipelotrichaceae were closely related to energy metabolism. Betaine correlated positively with [Eubacterium]_ventrio-sum_group, but negatively with Peptococcus and uncul-tured_bacterium_f_Erysipelotrichaceae, and a negative correlation was also detected between TMAO and un-cultured_bacterium _f_Erysipelotrichaceae. Therefore, the gut microbiota plays an important role in methylamine metabolism. TMAO is a metabolite derived from intestinal flora, while choline and betaine are the nutrient precursors of TMAO, they are each independently associated with the prevalence of CVD. Cholic acid and lithocholic acid shown strong negative to Treponema_2, but positive to Staphylococcus, Gordonibacter and Adlercreutzia. Cholic acid is one of the primary bile acids produced from cholesterol in the liver. Therefore, it can be considered H-HF is Fig. 6 Correlation analysis of relative abundance of gut microbiota in the genus level and fecal metabolite levels. Red means a positive correlation, and blue means a negative correlation. *P < 0.05, ** P < 0.01 related to dysbiosis of bile acid metabolism mediated by treponema_2, staphylococcus, gordonibacter and adlercreutzia.
These correlations indicated that perturbations in the gut microbiome, which may result in a significantly altered metabolomic profile; mainly including energy metabolism, amino acid metabolism, and methylamine metabolism. Our results further confirmed that the variations in gut microbiota and fecal metabolic phenotype associated with the development of HF.

Discussion
To the best of our knowledge, this is the first examination of specific changes in the gut microbiota composition and function in Dahl salt-sensitive rat model of hypertensive heart failure using both 16S rRNA and LC-MS metabolomics. Altered gut microbial composition and metabolites were observed in the H-HF rat model. The correlation analysis revealed that alterations in GM might contribute to H-HF through amino acid metabolism, bile acid metabolism, methylamine metabolism, energy metabolism and other aspects. The reduction in SCFA-producing bacteria and TMAO might be a notable characteristic for H-HF.
Recent evidence indicates that inflammation drives hypertensive heart failure [20], and impaired epithelial absorption may lead to translocation of microorganisms into the systemic circulation, possibly increasing HF by inducing systemic inflammation [21,22]. In our study, substantial inflammatory cellular infiltrate in the myocardial interstitium and integrity loss of the intestinal mucosal with infiltration of inflammatory cells into the colon tissue were observed in the H-HF rat model, which indicated inflammation involvement. Increased serum LPS and Zonulin indicated increased intestinal permeability and the loss of gut barrier function, which may lead to the development of metabolic endotoxemia and promote cardiac inflammation, thus exacerbating the development of heart failure.
Different from some clinical studies of GM [23,24], we found that alpha diversity analysis revealed no significant difference in gut microbiota diversity between groups based on Chao 1 and Shannon indices. However, another clinical study [11] showed the Chao 1 richness and the Shannon index were not significantly different between younger and older patients with HF, which is similar to our results.
To investigate taxonomic changes, we calculated the microbial abundance of the GM at the phylum level and family level. Altered gut microbiota composition was observed in H-HF group compared to CON and SR groups based on 16S rRNA gene sequencing results. The ratio of F/B increased at phylum level, while a decrease in the abundance of Muribaculaceae, Lachnospiraceae, and Lactobacillaceae was shown in the H-HF group compared with the CON group and SR group at family level, which suggested the occurrence of bacterial translocation in H-HF. Butyrate is an energy source for intestinal epithelial cells, modulating the epithelial barrier integrity and playing a local anti-inflammatory role in intestinal mucosa. Therefore, the decrease of butyrate content has been associated with an increase in endotoxemia and inflammation, and the decrease of butyrate-producing bacteria has been associated with inflammatory disease, including diabetes mellitus, obesity, hypertension, and inflammatory bowel disease [25]. Researchers have shown that the Lachnospiraceae family, which includes of several butyrateproducing species, is reduced in patients with heart failure [24]. The abundance of Muribaculaceae has also been proved to have a strong correlation with propionate and was an important predictor of SCFA concentrations [26]. As SCFA, the energy source for intestinal epithelial cells is also important in maintaining the integrity of intestinal epithelial cells, we could hypothesize that SCFA content is reduced due to the decrease of SCFA producing bacteria, and the imbalance of F/B ratio leads to abnormal energy [16,17],thus abnormal energy supply may directly affect the systolic function of the heart, which may lead to HF [18], the reduction in SCFA-producing bacteria might be a notable characteristic for H-HF, consistent with the results of clinical studies [27].
Moreover, Lactobacillus is essential in the maintenance of intestinal barrier function and integrity [28]. Studies have found that the decrease of Lactobacillus seems to promote the development of HF [29]. In our study, a decrease of Lactobacillus was also detected in H-HF groups. Lactobacilli can reduce cardiac hypertrophy and HF after myocardial infarction and improve left ventricular ejection fraction and shortening fraction, thus regulating intestinal flora may be used to develop a potential therapy to attenuate heart failure [30].
Comparing the different microbiotas between groups at the genus level, there were a total of 17 different flora among the three groups. The microflora can be considered as the biomarkers of hypertensive heart failure. To further assess the impact of a shifted gut microbiome on the host, we conducted metabolome analyses. A series of studies have been published which have highlighted the power of metabolic profiling to expand our understanding of metabolic derangements of the HF [31]. In our study, 35 metabolites have been identified as biomarkers for hypertensive heart failure, and 5 metabolic pathways were influenced, including: Histidine metabolism, Arginine and proline metabolism, Alanine, aspartate and glutamate metabolism, Glycine, serine and threonine metabolism and Glycerophospholipid metabolism. The pathways mainly involved energy metabolism and amino acid metabolism, it indicated that there are abnormities of energy metabolism and amino acid metabolism in H-HF.
Among these differential metabolites, TMAO is considered to be positively associated with cardiovascular disease. TMAO, a metabolites of the gut microbiota from specific dietary nutrients, is linked to a higher risk of HF, and a combination of TMAO and NT-proBNP could provide additional prognostic information [32]. Research shown that elevated plasma TMAO level in patients with HF is associated with poorer prognoses [33]. Systematic review and meta-analysis also demonstrated a positive dose-dependent relationship between TMAO plasma levels and increased cardiovascular risk and mortality [34]. Different from the present studies of TMAO, our study shown that the TMAO content in the feces of H-HF group was lower than the CON group and SR group, results were inconsistent perhaps owing to the different samples examined (e.g., plasma versus feces samples).
Gut microbiota perturbations associated with metabolic phenotype can be used to explore the possible mechanisms in the development of diseases. Peng et al., integrated 16S rRNA sequencing, metagenomics, and metabolomics to characterize gut microbial composition, function, and fecal metabolic phenotype in non-obese Type 2 diabetic goto-kakizaki rats [35], the results suggested that an altered gut microbiota is associated with T2DM pathogenesis. Yu et al. studied the variations in gut microbiota and fecal metabolic phenotype associated with depression by 16S rRNA gene sequencing and LC/ MS-based metabolomics [36], and showed a strong correlation between gut microbiota, fecal metabolites, and catecholamine levels. Herein, we observed a significant correlation between gut microbiota at genus level and fecal metabolites through Pearson's correlation analysis. The correlation indicated that perturbations in the gut microbiome, which may result in a significantly altered amino acid metabolism, bile acid metabolism, TMAO metabolism, energy metabolism and other aspects.
Aberrant energy metabolism is one of the most important signs of HF [37]. Dysfunction of energy metabolism can further promote the development and deterioration of HF. Therefore, an improvement in energy metabolism has been proposed as a potential treatment for HF. Recently, the importance of gut microbiota to the body's energy metabolism has been widely acknowledged. Our study shows that Adlercreutzia, Gordonibacter and uncultured_bacter-ium_f_Erysipelotrichaceae were closely related to energy metabolism. Therefore, regulation of related flora may be a new way to improve host energy metabolism.
Bile acids (BAs) are secreted by the liver and released into the small intestine during a meal where they aid in the absorption of dietary fat and fat-soluble vitamins [38]. BAs have emerged as important mediators of metabolic homeostasis, and have been proposed to play a direct role in regulating cardiovascular physiology and gut dysbiosis-induced HF [39]. BAs have been shown to regulate microbiota both directly and indirectly [40]. Our study shows that treponema_2, staphylococcus, gordonibacter and adlercreutzia were closely related to bile acid metabolism, which indicates that bile acid is closely related to intestinal flora, which is consistent with the literature.
Our results further confirmed that the variations in gut microbiota and fecal metabolic phenotype were associated with the development of HF, and an intervention to correct gut microbiota composition could be an innovative therapeutic strategy for HF.

Conclusion
In this study, an integrated approach of 16S rRNA gene sequencing combined with LC-MS based metabolomics was performed to assess the changes variations of gut microbiota in hypertensive heart failure rats. A total of 17 significantly altered bacterial genera and 35 metabolites were identified as the biomarkers of H-HF. Our results showed that HF significantly altered not only the gut microbiota composition but also fecal metabolic phenotype. In addition, correlation analysis revealed that some altered gut microbiota genera were strongly correlated with changed fecal metabolites. The reduction in SCFA-producing bacteria and TMAO might be a notable characteristic for H-HF. Overall, HF may contribute to changes in intestinal flora structure and metabolic function, regulated gut microbiota-related metabolites may be the potential biomarkers for diagnosis, prevention and treatment of HF.

Animals and treatments
Six-week-old Dahl salt-sensitive (SS, n = 16) rats and salt-resistant consomic SS.13BN (SR group, n = 8) rats with a body weight of 200-220 g were purchased from the Beijing Weitong Lihua Animal Co., Ltd., with a qualified number of 1,100,111,911,056,756,. Animals were housed in a specific pathogen-free area with ambient temperature (22 ± 5°C) and a 12/12 h light/dark cycle. After one week of adaptation, the SS rats were randomly allocated into two groups (8 per group): control (CON) group and hypertensive heart failure (H-HF) model. The CON group was given a low-salt diet containing 0.3% NaCl. To exclude the effects of high salt on the results, the SR group was set. The SR group and H-HF group were given a high-salt diet containing 8% NaCl for 20 weeks [41]. All feed was provided by Beijing Keao Xieli Feed Co., Ltd., (Beijing, China). Food and water were provided ad libitum throughout the experiments. Animal protocols in this study were supervised and approved by the Institutional Animal Care and Use Committee of Hunan University of Chinese Medicine.

Sample collections and preparation
The rats were sacrificed using urethane anesthesia (1.0 g/kg, i.p.). Fecal samples were collected immediately after defecation at the end of the experiment (20 weeks), and blood samples were collected from the abdominal aorta. All blood samples were processed into serum aliquots on the day of collection and stored at − 80°C just before the following analysis. Heart tissues were fixed in 4% paraformaldehyde solution, and 5-μm-thick paraffinembedded tissue sections were cut and stained with hematoxylin and eosin. The serum concentrations of NT-proBNP and LPS were detected using ELISA kits following the manufacturer's instructions (NT-proBNP Elisa Kit and LPS Elisa Kit: CUSABIO, Wuhan, China, Zonulin Elisa Kit: mlbio, Shanghai, China).

Echocardiography and blood pressure measurement
Cardiac functions were evaluated before the animals were sacrificed using an echocardiography method. After anesthetizing with urethane (1.0 g/kg, i.p.), all rats underwent echocardiography using a SonoScape-S2N ultrasound system (Shenzhen Kaili technology co., Ltd.). The following parameters were measured from twodimensional images and M-mode interrogation taken from the parasternal long-axis view at the papillary muscle level. The left ventricular ejection fraction (LVEF) and left ventricular fractional shortening (LVFS) were calculated according to the Teichholtz formula. Blood pressure was measured using a Volume Pressure Recording (VPR) system (CODA; Kent Scientific). For each animal, the systolic blood pressure (SBP) and diastolic blood pressure (DBP) were calculated as the average of 3 independent measurements.

Fecal metabolomics
Samples comprising 50 mg of fecal was placed in an EP tube, and 1 mL extraction solution (acetonitrile: methanol: water = 2:2:1, with isotopically-labelled internal standard mixture) was added. After 30s vortexing, the samples were homogenized at 35 Hz for 4 min and sonicated for 5 min in ice-water bath. The homogenization and sonication cycle were repeated for 3 times. Then the samples were incubated for 1 h at − 40°C and centrifuged at 12000 rpm for 15 min at 4°C. The resulting supernatant was transferred to a fresh glass vial for analysis. The quality control (QC) sample was prepared by mixing an equal aliquot of the supernatants from each samples (Fig. S1a). LC-MS/MS analyses were performed using an UHPLC system (Vanquish, Thermo Fisher Scientific) with a UPLC BEH Amide column (2.1 mm × 100 mm, 1.7 μm) coupled to Q Exactive HFX mass spectrometer (Orbitrap MS, Thermo). The raw data were converted to the mzXML format using ProteoWizard and processed with an in-house program, which was developed using R and based on XCMS, for peak detection, extraction, alignment, and integration. Then an inhouse MS2 database (BiotreeDB) was applied for metabolite annotation. The cutoff for annotation was set at 0.3. The methods used in this study were in accordance with the published literature [13,42]. Data was acquired in positive and negative ion modes, with the two sets of data combined for analysis.
In this study, 13,929 peaks in positive ion mode and 11,910 peaks in negative ion mode were detected, and among them, 1066 metabolites were found in positive ion mode, while 346 metabolites were found in negative ion mode after relative standard deviation de-noising.
The data were trimmed using Compound Discoverer 2.1 (Thermo Fisher Scientific, Waltham, MA, United States) and imported into SIMCA16.0.2 software package (Sartorius Stedim Data Analytics AB, Umea, Sweden) for principle component analysis (PCA) and orthogonal projections to latent structures discriminate analysis (OPLS-DA). Then, a 7-fold cross validation was performed to calculate the value of R 2 and Q 2 .Furthermore, the value of variable importance in the projection (VIP) of the first principal component in OPLS-DA analysis was obtained. It summarizes the contribution of each variable to the model. The metabolites with VIP > 1.0 and p < 0.05 (Student's t test) were considered as significantly changed. In addition, commercial databases including KEGG (http://www.genome.jp/kegg/) and MetaboAnalyst (http://www.metaboanalyst.ca/) were used for pathway enrichment analysis. The Data analysis method used in this study is consistent with published literature [42].

16S rRNA gene sequencing analysis
Total genomic DNA from fecal samples was extracted by Tiangen Fecal Genomic DNA Extraction Kit under the manufacturer's instruction. The V3-V4 regions of 16S rRNA genes were PCR-amplified using the following primers: 338F: 5′-ACTCCTACGGGAGGCAGCA-3′ and 806R: 5′-GGACTACHVGGGTWTCTAAT-3′, and amplification products were purified, quantified and homogenized to get a sequencing library. Library QC was performed for constructing libraries, qualified libraries were sequenced on Illumina HiSeq 2500. The original image data files obtained by Illumina HiSeq highthroughput sequencing were converted into Sequenced Reads by Base Calling analysis, the results were stored in FASTQ format files. Paired-ends sequences were merged then filtered in three steps: 1) PE reads merge: FLASH v1.2.7 [43] software was used to merge reads through overlap, the obtained merged sequences were Raw Tags; 2) Tags filtering: Trimmomatic v0.33 software was used to filter merged Raw Tags to get high quality Clean Tags;3) Remove Chimera: UCHIME v4.2 software was used [44] to identify and remove chimeric sequences to get Effective Tags.
Sequence analysis was performed by Uparse software (Uparse v7.0.100, http://drive5.com/uparse/) [45]. Sequences with ≥97% similarity were assigned to the same OTU. Representative sequence for each OTU was compared to the Silva Database (http://www.arb-silva.de/) [46] using Mothur (version v.1.30) to identify taxonomic information. In order to study phylogenetic relationships between OTUs, multiple sequence alignment was conducted using the MUSCLE software (Version 3.8.31, http://www.drive5.com/muscle/) [47] . All samples were normalized and mothur software was used to analyze alpha diversity of samples. QIIME software (Version 1.9.1) was used to perform beta diversity analysis. Finally, ANOVA analysis was performed to screen species with significant differences at genus level (P < 0.05).

Statistical analysis
Statistical analysis was performed using IBM SPSS Statistics 22.0 (Chicago, USA). Differences between groups were evaluated by one-way analysis of variance (ANOVA). The relationships between fecal metabolites and the gut microbiome taxa were assessed using Pearson correlation. The significance threshold was set at P < 0.05 for all tests.