Exploration of the Characteristics of Intestinal Microbiota and Metabolomics in Different Rat Models of Mongolian Medicine

Background Mongolian medicine is a systematic theoretical system, which is based on the balance among Heyi, Xila, and Badagan. However, the underlying mechanisms remain unclear. This study aimed to explore the characteristics of intestinal microbiota and metabolites in different rat models of Mongolian medicine. Methods After establishing rat models of Heyi, Xila, and Badagan, we integrated 16S rRNA gene sequencing and metabolomics. Results Heyi, Xila, and Badagan rats had significantly altered intestinal microbial composition compared with rats in the MCK group. They showed 11, 18, and 8 significantly differential bacterial biomarkers and 22, 11, and 15 differential metabolites, respectively. The glucosinolate biosynthesis pathway was enriched only in Heyi rats; the biosynthesis of phenylpropanoids pathway and phenylpropanoid biosynthesis pathway were enriched only in Xila rats; the isoflavonoid biosynthesis pathway, the glycine, serine, and threonine metabolism pathway, and the arginine and proline metabolism pathway were enriched only in Badagan rats. Conclusions The intestinal microbiota, metabolites, and metabolic pathways significantly differed among Heyi, Xila, and Badagan rats compared with control group rats.


Introduction
Traditional Mongolian medicine is an indigenous medicine system widely practiced in China, especially in the Inner Mongolia region [1]. e systematic theoretical system of Mongolian medicine is based on the balance among three roots: Heyi, Xila, and Badagan. Generally, the ratio of three roots in different individuals depends on genetic and environmental factors. An imbalance in these roots results in disease. e Heyi, Xila, and Badagan rat models constructed based on the "Four-Part Medicine Classics" [2] would enhance our understanding of Mongolian medicine. Current research on Mongolian medicine has mainly focused on clinical practice or drugs prescribed in Mongolian medicine [3][4][5]. However, the underlying mechanisms based on Heyi, Xila, and Badagan are unclear. e digestive system status usually plays an important role in the diagnosis of different diseases in Mongolian medicine, among which the gut environment and intestinal microbiota carry significant importance [6,7]. erefore, it is important to understand the underlying pathogenesis of different diseases by comparing the composition of intestinal microbiota in Heyi, Xila, and Badagan rat models. e intestinal microbiota refers to the various microorganisms present in the gastrointestinal tract, including bacteria, fungi, and viruses [7]. About 3 × 10 6 genes exist in microbial genome sequences [8]. With the development of next-generation sequencing technology, 16S rRNA sequencing made it possible for us to link intestinal microbiota to various diseases [9]. Dysbiosis is known to cause diseases, such as hypertension [10], inflammatory bowel diseases [11], and type 2 diabetes [12]. To the best of our knowledge, no previous study has linked intestinal microbiota to the different aspects of Heyi, Xila, and Badagan in Mongolian medicine. After oral administration, drugs used in traditional medicine often work by interacting with the intestinal microbiota [9,13], implying that the intestinal microbiota is important in traditional Chinese medicine and Mongolian medicine. Metabolomics analysis has emerged as an effective tool to study pathogenetic mechanisms [14]. Some recent studies have combined metabolomics analysis with intestinal microbiota analysis [15][16][17]. A recent study in a twin model demonstrated complex links between host phenotypes and intestinal microbiota based on metabolic profiling [18].
us, metabolomics analysis and determination of the composition of intestinal microbiota would help us better understand different diseases in Mongolian medicine.
In the current study, we explored the characteristics of intestinal microbiota and metabolites via an integrated analysis of 16S rRNA sequencing and metabolomics in three Mongolian medicine rat models (Heyi, Xila, and Badagan rat models) to further understand the possible mechanisms underlying diseases related to the three roots in Mongolian medicine.
e results would also provide deeper insights into the function of intestinal microbiota. ere was no intervention put on rats in the control group (MCK). And all rat models were constructed mainly based on the "Four-Part Medicine Classics" [2]. When the rats showed corresponding characteristics of Heyi, Xila, and Badagan described in the "Four-Part Medicine Classics," the model was considered to be constructed successfully (Supplementary Material 1). A Heyi rat model was constructed according to the following methods: First, diet intervention, drinking water was replaced by cold black tea (5 g tea + 100 mL distilled water) and the rats were fed buckwheat (8.5 g/day). Second, behavior intervention, the rats were exposed to the continuous cat audio at 70 decibels. ird, Mongolian medicine intervention, the rats were given a dose of 1 mL/100 (g d) Gaburi by gavage. Finally, 0.1 mL tail vein bloodletting was performed on the rats at 5 pm every two days. It took 31 days to construct the Heyi rat model. A Xila rat model was constructed according to the following methods: First, diet intervention, the rats were given 1 mL liqueur by gavage once every other day, were given 0.7 g/kg fruit oil at 6 am every day, and were fed yellow rice (15 g/day). Second, behavior intervention, the rats were under the environment of 29 ± 2°C. ird, Mongolian medicine intervention, the rats were given 0.7 g/kg pepper by gavage daily at 12 noon. It took 21 days to construct the Xila rat model. A Badagan rat model was constructed according to the following methods: First, diet intervention, the rats were fed lard and wheat flour (mixed in a ratio of 1 : 4). Second, behavior intervention, the rats were under the environment of 60 ± 5% humidity.

Materials and Methods
ird, Mongolian medicine intervention, the rats were given 4 mL dandelion (200% decoction) by gavage. It took 49 days to construct the Badagan rat model.

Fecal Sample Collection and DNA Extraction.
Fecal samples were collected in a sterile conical tube and stored at −80°C. According to the manufacturer's instruction, the DNA was extracted using an E.Z.N.A. feces DNA kit (Omega Bio-Tek, Norcross, GA, USA). e DNA quality was determined by 1% agarose gel electrophoresis, and the DNA concentration was determined using a NanoDrop 2000 spectrophotometer ( ermo Fisher Scientific, USA).

16S rRNA Microbial Community Analysis.
e primer 341F (5′-CCTAYGGGRBGCASCAG-3′) 806R (5′-GGAC-TACNNGGGTATCTAAT-3′) was used to amplify the bacterial 16S rRNA gene V3-V4 region, which was performed on the Illumina HiSeq sequencing platform (Illumina, USA). As for the paired-end sequences obtained, the primer adapter sequences were removed, and various samples were distinguished based on the barcode tag sequences. And the valid data of the samples were obtained after the quality control filtering. FLASH software (version 1.2.11) [19] was used to splice the paired-end sequences, and Trimmomatic software (version 0.33) [20] was used to filter the spliced sequences. UCHIME software (version 8.1) [21] was used to remove the chimera sequences in order to obtain valid data for further analysis. Based on 97% similarity, all sequences were clustered in operational taxonomic units (OTU) using USEARCH software (version 10.0) [22], which was filtered with 0.005% of all sequences as a threshold. In order to determine the classification, RDP Classifier software (version 2.2) (http://rdp.cme.msu.edu/classifier/classifier. jsp) [23] was used to compare the representative sequence of each OTU with the Silva database (https://www.arb-silva. de/) [24]. e alpha diversity of microbiota was calculated using mothur software (version 1.30) [25], including Chao1, ACE, Shannon index, and Simpson index. e β diversity was estimated according to the Bray Curtis distance algorithm and then was visualized using nonmetric multidimensional scaling (NMDS). e differential biomarkers between different groups were found according to the linear discriminant analysis effect size (LEfSe) [

Serum Sample Preparation for Metabolome.
Blood samples were collected from the abdominal aorta, and the serum was separated and stored at −80°C. en, 200 μL of the serum was taken and 3 times volume of precooled acetonitrile solution was added. e sample was vortexed and mixed, and then it was placed in a refrigerator at −20°C for 30 min. Subsequently, the sample was centrifuged (14000 g, 4°C, 15 min) and the supernatant was transferred in a new centrifuge tube for concentration and drain. e 1 : 1 (v/v) mixture of mobile phase A (ammonium acetate) and mobile phase B (acetonitrile) was used to redissolve the sample; after the high-speed centrifugation, the sample was used for HPLC-MS analysis. e target compounds were separated on an Accucore Hilic (100 × 2.1 mm, 2.6 μm) liquid chromatography column, using Vanquish ( ermo Fisher Scientific) ultra-performance liquid chromatography. And the liquid chromatography consisted of 10 mM ammonium acetate as phase A and acetonitrile/10 mM ammonium acetate as phase B (9 : 1). Gradient elution was used: 0∼1 min, 100% A; 1∼9 min, 0%∼100% B; 9∼12 min, 100% B; and 12.1∼15 min, 100% A. e flow rate of mobile phase was 0.35 mL/min; the column temperature was 35°C; the sample tray temperature was 4°C; and the injection volume was 2 μL.

Metabolite Analysis.
e serum metabolites were analyzed using ermo Scientific's Q Exactive mass spectrometer. e positive and negative ions were scanned once each. e positive ion scan was performed firstly, after which, the negative ion scan was performed. e full scan range was 70-1050 m/z. In the full scan, the precursor ions with TOP10 ion intensity were selected for secondary MS identification. e precursor ion was fragmented according to the HCD method, which was used for secondary mass spectrometry sequence determination and then generated the mass spectrometry detection original file. Subsequently, the raw data were transformed in mzML format using ProteoWizard software, and XCMS was used to perform retention time correction, peak identification, peak extraction, peak integration, peak alignment, etc. en, identification of metabolites was based on the Compound Discover V3.0 (CD) and mzCloud database. SIMCA-P software was used for orthogonal partial least square discriminate analysis (OPLS-DA). In order to select different variables as potential markers, VIP-plot (VIP >1) was obtained from OPLS analysis. e differential metabolites with VIP >1 and P value < 0.05 were screened.
e OPLS-DA model was validated using 7-fold crossvalidation. en, R 2 Y (model explainability of the categorical variable Y) and Q 2 (predictability of model) were used to determine the validity of the model. Finally, the permutation test was performed to further test the validity of the model, which was done via randomly changing the arrangement of categorical variable Y (n � 200 times) and obtaining random Q 2 values.

Statistical Analysis.
e Wilcoxon rank-sum test (R software v3.6.2) was used to compare the α diversity (ACE index, Chao1 index, Shannon index, and Simpson index) and microbiota between various groups, and P < 0.05 was considered as the significance threshold. Analysis of similarities (ANOSIM) was used to analyze the differences between and within groups. e Kruskal-Wallis sum-rank test was used to determine the alterations in abundance between different groups in LEfSe analysis, and |LDA score| > 3 and P < 0.05 were taken as the difference screening thresholds.

Changes in Intestinal Microbiota
Diversity among Different Mongolian Medicine Rat Models. Based on the results of 16S rRNA sequence analysis, the changes of intestinal microbiota of three different Mongolian medicine rat models and control group rats were investigated. e results were clustered in operational taxonomic units (OTU) based on over 97% similarity. e rarefaction curves, based on the number of sample reads and OTUs, tended to be flat ( Figure S1), indicating that the amount of sequencing data was sufficient to reflect the species diversity in all samples. e ACE, Chao1, Shannon, and Simpson indexes were used to evaluate microbial alpha diversity.
Compared with the MCK group, the ACE and Chao1 indexes of intestinal microbiota in the Heyi rat model significantly decreased (P value < 0.05), but the Shannon and Simpson indexes did not significantly differ. Compared with the MCK group, the ACE and Chao1 indexes of intestinal microbiota in the Xila rat model did not significantly differ, but the Shannon index increased significantly and the Simpson index decreased significantly. Compared with the MCK group, the ACE index of intestinal microbiota in the Badagan rat model significantly decreased, but the Chao1, Shannon, and Simpson indexes did not significantly differ (Figures 1(a)-1(d)). e results showed that compared with the MCK group, the richness of intestinal microbiota in Heyi rats and Badagan rats was decreased, but the diversity did not significantly differ; however, the richness of intestinal microbiota in Xila rats did not significantly differ, but the diversity increased.
Compared with Heyi and Badagan rat models, the ACE, Chao1, and Shannon indexes of intestinal microbiota in the Xila rat model significantly increased, but the Simpson index decreased significantly. Compared with the Heyi rat model, the Shannon index of intestinal microbiota in the Badagan rat model decreased significantly, but the ACE, Chao1, and Simpson indexes did not significantly differ (Figures 1(a)-1(d)).
e results above indicated that the richness and diversity of intestinal microbiota in Xila rats both increased significantly compared with Heyi and Badagan rats. Moreover, compared with Heyi rats, the richness of intestinal microbiota in Badagan rats did not significantly differ, but the diversity decreased significantly.
According to the results of β diversity analysis, significant differences were noted between rats in the Heyi and MCK groups ( significantly different from that of the MCK group. Collectively, the intestinal microbial structure of Heyi, Xila, and Badagan rats significantly differed from that of the control group rats.
To investigate the specific intestinal bacterial biomarkers at the genus level, line discriminant analysis (LDA) effect size (LEfSe) analysis was performed on all three different Mongolian medicine rat models. e bacteria abundance of 23 genera in the Heyi rat model was significantly higher than that in MCK rats (LDA >3, P value < 0.05) (Figure 3(d)). e bacteria abundance of 30 genera in the Xila rat model was significantly higher than that in MCK rats (LDA >3, P value < 0.05) (Figure 3 index. X axis: different groups; Y axis: the corresponding diversity index. e Wilcoxon test was used to determine the significance of differences between any two groups. Statistical significance: * P < 0.05, * * P < 0.01, * * * P < 0.001, and * * * * P < 0.0001. genera in the Badagan rat model was significantly higher than that in MCK rats (LDA >3, P value < 0.05) (Figure 3( (Figure 3(g)). Compared with the MCK group, the abundance of 3 KEGG pathways significantly increased and that of 10 pathways significantly decreased in the Heyi rat model, total 13 KEGG pathways were significantly different between Heyi and MCK groups (Figure 4(a)). Compared with the MCK group, the abundance of 4 KEGG pathways significantly increased and that of 13 pathways significantly decreased in the Xila rat model, and 17 KEGG pathways were significantly different between Xila and MCK groups (Figure 4(b)). Compared with the MCK group, the abundance of 6 KEGG pathways significantly increased and that of 14 pathways decreased significantly in the Badagan rat model, and 20 KEGG pathways were significantly different between the Badagan and MCK groups (Figure 4(c)).

Prediction of the Function of Intestinal Microbiota in
Moreover, the abundance of 4 KEGG pathways was significantly different between all three experimental groups (Heyi, Xila, and Badagan rat models) and the control group Y axis: beta distance. All between: beta distance of samples in all groups; all within: beta distance of samples within the group. P value < 0.05 was considered as significant difference. R > 0 indicated that the difference between groups was greater than the difference within the groups; R < 0 indicated that the difference within the group is greater than the difference between the groups; the greater the |R| value, the greater the relative difference. ( Figure 4(d)), among which the abundance of "global and overview maps" pathway and "energy metabolism" pathway was significantly increased in the experimental groups ( Figure 4(e)), but that of "digestive system" pathway and "endocrine and metabolic diseases" pathway was significantly reduced (Figure 4(f )).     Furthermore, the abundance of some pathways was only changed in one certain model. e abundance of "cardiovascular diseases" pathway and "neurodegenerative diseases" pathway was significantly decreased only in the Heyi rat model. e abundance of "metabolism of cofactors and vitamins" pathway was significantly increased and the abundance of "nucleotide metabolism" pathway, "replication and repair" pathway, "infectious diseases: parasitic" pathway, and "cancers: overview translation" pathway was decreased only in the Xila rat model. e abundance of "cellular communityprokaryotes" pathway was significantly increased and the abundance of "transport and catabolism" pathway, "glycan biosynthesis and metabolism" pathway, "cell motility" pathway, "environmental adaptation" pathway, "excretory system" pathway, "aging" pathway, "immune system" pathway, "signal transduction" pathway, and "membrane transport" pathway was decreased only in the Badagan rat model. Our data implied that the abundance of various pathways has been altered in different models compared with the MCK group.

Identification of Serum Metabolic Profile and Metabolic Markers in Different Mongolian Medicine Rat Models.
According to the results of metabolite profile analysis, Heyi, Xila, and Badagan rat models presented significant differences compared with the MCK group (Figures 5(a)-5(f )).
In the OPLS-DA multivariate model, metabolites with VIP scores >1 and P value < 0.05 were considered as differential metabolites. Compared with the MCK group, 30 differential metabolites were detected in the Heyi rat model after positive ion scan and 64 differential metabolites were noted after negative ion scan (Table S4), resulting in 94 differential metabolites involved in 7 metabolic pathways ( Figure 5(g)). Compared with the MCK group, 35 differential metabolites were detected in the Xila rat model after positive ion scan and 51 differential metabolites were noted after negative ion scan (Table S5), resulting in 86 differential metabolites involved in 8 metabolic pathways ( Figure 5(h)). Compared with the MCK group, 37 differential metabolites were detected in the Badagan rat model after positive ion scan and 63 differential metabolites were noted after negative ion scan (Table S6), resulting in 100 differential metabolites involved in 8 metabolic pathways (Figure 5(i)).
In addition, 22, 11, and 15 differential metabolites were detected only in the Heyi, Xila, and Badagan rat models, respectively ( Figure 5(j), Table S7). Moreover, the "glucosinolate biosynthesis" pathway was enriched only in the Heyi rat model; the "biosynthesis of phenylpropanoids" pathway and "phenylpropanoid biosynthesis" pathway were enriched only in the Xila rat model; the "isoflavonoid biosynthesis" pathway, "glycine, serine, and threonine metabolism" pathway, and "arginine and proline metabolism" pathway were enriched only in the Badagan rat model.

Discussion
We explored the characteristics of intestinal microbiota and metabolomics in different Mongolian medicine rat models. In this study, a joint analysis of 16S rRNA gene sequencing and metabolomics was conducted to investigate the potential mechanisms underlying diseases related to Heyi, Xila, and Badagan. Our results indicated that the intestinal microbiota of Heyi, Xila, and Badagan rat models significantly differed from that of control rats. Metabolites and metabolic pathways in Heyi, Xila, and Badagan rats were also significantly different from those in control rats. e alpha and beta diversity and intestinal microbiota composition were investigated, as the link between intestinal microbiota changes and many metabolic diseases have been reported [28,29]. Our results showed that compared with the MCK group, the alpha diversity of intestinal microbiota in Xila rats was increased, but that in Heyi and Badagan rats showed no significant difference. e alpha diversity of intestinal microbiota in Xila rats was also higher than that in Heyi and Badagan rats. e results of beta diversity indicated that there was a significant dissimilarity between the control group and Mongolian medicine rat models. e intestinal microbiota of Heyi, Xila, and Badagan rats was significantly altered compared with that of control group rats, which was probably an important contributor to different diseases in Mongolian medicine. Furthermore, intestinal microbiota composition was also investigated. Firmicutes and Bacteroidetes were found to be the dominant strains in all rats, consistent with former studies reporting that Firmicutes and Bacteroidetes are two main phyla comprising gut microbiota in healthy humans [30]. Moreover, the abundance of some phyla in our disease models was significantly different from that in control group rats. Verrucomicrobia, usually colonized in the mucosal layer, is considered a promising probiotics [31,32]. e abundance of Verrucomicrobia was increased in Heyi rats but decreased in Badagan rats compared with control group rats, which might be responsible for the contrary manifestations of Heyi disease and Badagan disease. e abundance of Proteobacteria was increased both in Xila and Badagan rats, and it is reported that the increase in abundance of Proteobacteria resulted in an imbalanced gut microbiota composition and consequently metabolic disorders [33]. Furthermore, we found that bacteria biomarkers of some genera existed specifically in certain disease in Mongolian medicine. ere were 11, 18, and 8 bacterial biomarkers that were increased only in Heyi, Xila, and Badagan rats, respectively. Some of these bacterial genera have been associated with diseases, such as Faecalibacterium [34], Prevotella_1 [35], Prevotellaceae [36], Alistipes [37], Bilophila [38], Hungatella [39], and Sellimonas [40], which might contribute to diseases in Mongolian medicine. In addition, the differential KEGG pathways were found between Mongolian medicine rat models and control group rats. e abundance of "global and overview maps" pathway and "energy metabolism" pathway was significantly increased in the disease models, but the abundance of "digestive system" pathway and "endocrine and metabolic diseases" pathway was significantly decreased. We suspected that these pathways might play an essential role in three Mongolian medicine rat models. However, further studies should be conducted to further understand the intestinal microbiota in different diseases in Mongolian medicine.

Evidence-Based Complementary and Alternative Medicine
Furthermore, investigations of metabolites and metabolic pathways in Heyi, Xila, and Badagan rats revealed 22, 11, and 15 differential metabolites specific to Heyi, Xila, and Badagan rats, respectively. e "glucosinolate biosynthesis" pathway was enriched only in Heyi rats, whereas the "biosynthesis of phenylpropanoids" pathway and "phenylpropanoid biosynthesis" pathway were enriched only in Xila rats.
e "isoflavonoid biosynthesis" pathway, "glycine, serine, and threonine metabolism" pathway, and "arginine and proline metabolism" pathway were enriched only in Badagan rats. ese specific metabolic pathways probably play an important role in the characteristics of the various Mongolian medicine rat models. Moreover, we noticed that most of these pathways were related to the metabolism of certain amino acids, which might be correlated with the diseases in Mongolian medicine. For example, it has been reported that glycine-conjugated metabolites were involved in chronic kidney disease and hypertension in rats [41]. Gut-derived D-serine has renoprotective effects on the kidney in acute kidney injury [42]. L-arginine protects the intestinal barrier by promoting expression of tight junction proteins in rats [43]. erefore, the kidney and gut are probably more important for the manifestations of diseases in Mongolian medicine. Collectively, although metabolomics analysis helps us better understand the three disease aspects in Mongolian medicine, the specific role of each pathway in various diseases in Mongolian medicine remains to be further studied.
In spite of this, there are still several limitations in our present study. First, only 16S rRNA sequencing and metabolomics were included in our research, which might lead to inevitable deviation in the results. Multiple omics data like methylation data could be further studied in the future. Moreover, the detailed underlying reasons for the pathway changes should be further explored.

Conclusions
In conclusion, we have firstly explored the characteristics of intestinal microbiota and metabolomics in different Mongolian medicine rat models by integrating 16S rRNA sequencing and metabolomics approaches. Our data showed that the intestinal microbiota, metabolites, and metabolic pathways in Heyi, Xila, and Badagan rats were significantly different from those in control group rats. Despite the detailed mechanisms remain to be clarified, our research has provided more reference information for diseases in Mongolian medicine.

Data Availability
e data sets of this study are available on request to the corresponding author. is study was also supported by the National Natural Science Foundation of China (No. 82060912, to Eerdunchaolu). Figure S1: the rarefaction curves of all samples. Table S1: relative abundance of microbial phylum (percentage) in the Heyi rats and control rats. Table S2: relative abundance of microbial phylum (percentage) in the Xila rats and control rats. Table S3: relative abundance of microbial phylum (percentage) in the Badagan rats and control rats. Table S4: differential metabolites of Heyi rat samples compared with control group. Table S5: differential metabolites of Xila rat samples compared with control group. Table S6: differential metabolites of Badagan rat samples compared with control group. Table S7: differential metabolites only present in a group of rats. (Supplementary Materials)