Tibetan wild asses (Equus Kiang) were superior to domestic donkeys (Equus Asinus) in terms of the composition and function of gut microbial community

Background: Tibetan wild asses are the only wild species of perissodactyls on the Qinghai-Tibetan Plateau, and appears on the International Union for Conversation of Nature (IUCN) 2012 Red List of threatened species. The gut microbiota has a great effect on the health and nutrition of the host, however, scant research is available on the characteristics of their intestinal microbiota. Therefor, understanding the gut microbita composition and function of TWAs can provide a theoretical for the situ conservation of wild animals in the future. Results: To characterize its composition and function, we analyzed the intestinal microbiota of wild asses and domestic donkeys by high-throughput sequencing of the 16s rDNA regions. No significant difference in alpha diversity was detected between these two groups. Beta diversity showed that the bacterial community structure of wild asses was acutely different from domestic donkeys. At the phylum level, the two dominant phyla of Bacteroidetes and Firmcutes in wild asses were significantly higher than that in domestic donkeys. At the genus level, Ruminococcaceae_NK4A214 , Phascolarctobacterium , Coprostanoligenes_group , Lachnospiraceae_XPB1014_group and Akkermansia in wild asses were significantly higher than domestic donkeys. Moreover, statistical comparisons showed that 40 different metabolic pathways exhibited significant differences . Among them, 29 pathways had richer concentrations in wild asses than domestic donkeys, mainly amino acid metabolism, carbohydrate metabolism, and energy metabolism. Of note, network analysis showed that wild asses harbored a relatively more complex bacterial network than domestic donkeys, possibly reflecting the specific niche adaption of gut bacterial communities through species interactions. Conclusions: Wild asses were superior to that of domestic in gut microbial community composition and function. For wild animal conservation, wild asses are more suitable to wild than to be

herbivores, donkeys utilize forage nutrition by microorganism fermentation in the gastrointestinal tract. Zhao et al. (2016) studied bacterial diversity in the ceca of Xinjiang donkeys and found that Ruminococcaceae and Lachnospiraceae played a key role in digesting roughage feed [11]. Due to the difficulty of obtaining samples from wild protected animals, research on gut microbes in TWAs is limited. A recent study compared the gut microbiota composition of wild and captive TWAs, revealing that captivity reduced intestinal microbial diversity and increased the risk of epidemics, thereby negatively impacting the health of wild life [12]. However, few studies have compared the intestinal microbial composition and function of TWAs and natural pasture domestic donkeys (NPDDs) on the QTP.
As members of the same genus-Equus-,, TWAs and NPDDs provide ideal candidates for the study of gut microbiota. Therefore, we used high-throughput sequencing of the V3-V4 regions of the 16s rRNA gene to investigate the composition and predict the metabolic pathways and functions of the gut microbes in TWAs and NPDDs. We hypothesized that TWAs were superior to NPDDs in terms of the composition and function of their intestinal flora. We think our research is of great significance to wildlife conservation. 12 923 operational taxonomic units (OTUs) at a cutoff of 97% similarity were generated.
Among these OTUs, the number of unique OTUs in TWAs and NPDDs was 2948 and 2816, respectively, and 7159 OTUs were shared by all samples, accounting for 55.40% of the total OTUs.

Diversity analysis of gut microbiota in TWAs and NPDDs
To determine alpha diversity, we calculated the Chao1, ACE, observed species estimates and the Shannon diversity index. As shown in Fig. 1,no significant differences were found in these two groups for the Shannon index, Chao1, ACE estimates and the number of observed species. The OTU level rarefaction curves of diversity estimators reached a plateau at the minimum sequence number of 50,000 tags (Fig. S2).. With regard to beta diversity, PCoA analysis showed that the two groups of bacterial communities differed significantly ( Fig. S3A) ( P < 0.05). To assess the overall similarity, UPGAMA with classification tree analysis was conducted as shown in Fig. S3B. It is obvious that the bacterial community of TWAs was completely different from that of NPDDs.

Significance analysis of gut bacterial community abundance between TWAs and NPDDs
At the phylum level illustrated in Fig. 2A, Bacteroidetes, Firmicutes, Verrucomicrobia, and Fibrobacteres constituted the four dominant phyla in all samples. The average ratio of Firmicutes/Bacteroidetes was 0.80 and 0.74, respectively. The significance analysis of the top 10 phyla of microbial communities between TWAs and NPDDs is shown in Fig. 2B. The relative abundance of Bacteroidetes, Firmicutes, Cyanobacteria, and Synergistetes in TWDs was significantly higher than in NPDDs, whereas Lentisphaerae in TWAs was significantly lower (P < 0.05). No significant differences were detected for the relative abundances of Verrucomicrobia, Fibrobacteres, Spirochaetae, Proteobacteria, and detected in the gut samples (Table S2).. As shown in Fig. 4B, 18 gene families were remarkably different and the relative abundance of amino acid metabolism, energy metabolism, glycan biosynthesis, and the biosynthesis of other secondary metabolites in TWAs were significantly higher than in NPDDs (P < 0.05). As for membrane transport and cell motility, NPDDs had higher expression abundances than TWAs (P < 0.05). At KEGG level 3, the principal component analysis and partial least squares discriminant analysis revealed clear differences between the two groups (Fig. 5). As shown in Table S3, 40 metabolite pathways were significantly different (VIP > 1 and P < 0.05), of which ten metabolites pathways belonging to amino acid metabolism were enriched in TWAs (P < 0.05). A total of ten metabolic pathways belonged to carbohydrate metabolism of which seven metaboliteswere enriched in TWAs, whereas the remaining three metabolites (starch and sucrose metabolism, amino sugar and nucleotide sugar metabolism, and lipopolysaccharide biosynthesis) were enriched in NPDDs (P < 0.05). In addition, five metabolic pathways belonged to metabolism of cofactors and vitamins, of which two metabolites were higher in TWAs. For the remaining 15 metabolic pathways, 11 had a higher relative abundance in TWAs.

Network analysis of bacterial communities and putative keystone taxa
The putative keystone genus in bacterial networks were identified based on the Pi and Zi values (Table 2). Three key module hubs and two connectors were identified in NPDD.
Among these OTUs, module hubs were affiliated with Firmcutes, Cyanobacteria and unclassified bacteria, and the connectors were Anaerovibrio and unclassified bacteria. In contrast, four connector and no module hubs were identified in TWAs. Among these OTUs, two connectors were affiliated with Firmcutes and Actinobacteria and the others were the genus Oscillospira and Arthrobacter. Notably, these keystone genus in network showed relatively low abundances in gut bacterial communities.

Discussion
Although TWAs and NPDDs are herbivores species on the QTP, scant research has been conducted on their gut microbiota diversity and composition. To the best of our knowledge, this is the first study to compare the gut microbiota composition and function of TWAs and NPDDs. We found that the Shannon-Wiener index and ACE and Chao1 estimates, as well as the number of observed species between TWAs and NPDDs were not significantly different, suggesting that altitude and grassland types did not affect the alpha diversity of gut microbes. A previous study had shown that diet was one of the factors that affected gut microorganism diversity [13]. In our study, there was no significant difference in herbage nutrient composition between TWAs and NPDDs in different meadow types (swamp meadow and alpine meadow, respectively), which was consistent with the results of the alpha diversity analysis. In terms of biological classification, TWAs and NPDDs are animals of the same genus but different species; however, there was no significant difference in the alpha diversity of gut microbes, suggesting that species differences did not affect the alpha diversity of bacteria community, as was also found in a previous study on wild and domestic yak [14]. As to beta diversity, PCoA results indicated that the composition of gut microbiota was significantly different in these two groups. For a better understanding of the similarity of their gut microbiota, UPGMA methods were applied and the diagram revealed that the fecal microorganism communities in TWAs occupied their own phylogenetic tree branch which clearly distinguished from the fecal microorganisms in NPDDs. Therefore, we concluded that species, altitude, and meadow type can affect the intestinal bacterial structure; however, the determination of which of these played the leading role in shaping intestinal flora requires further studies.
To assess the composition of bacteria community, we first analyzed at the phylum level. In our study, Bacteroidetes was the most abundant phylum followed by Firmicutes in both TWAs and NPDDs, constituting more than 80% of the total bacteria content. This has been confirmed by former studies on intestinal microbial diversity in mammals [15][16][17], as these microorganisms facilitate the digestion of cellulose and hemicellulose in the diet [18]. The relative abundance of Bacteroidetes and Firmicutes in TWAs was significantly higher than in NPDDs, indicating that the cellulose and hemicellulose decomposition capacity was stronger in TWAs than in NPDDs under similar forage nutritional composition.
Moreover, the ratio of Bacteroidetes and Firmicutes affected the energy acquisition and body fat accumulation in humans and mice [19]. Previous studies have shown that a decrease in Bacteroidetes was strongly linked to an increase of fat in host tissue in mice and humans [19,20]. However, scant research has focused on the relationship between gut microbiota and obesity in TWAs and NPDDs. By observing the body size of TWAs and NPDDs, we speculated that the lean body size of TWAs might be related to the high Bacteroidetes content in their intestinal microbiota. At the genus level, 15 core genera whose relative abundances were more than 1% were detected in all gut samples; these genera constituted 32.43% of the total genera reads and were mainly affiliated with the phylum Firmicutes (nine genera) and phylum Bacteroidetes (four genera). The others were affiliated with the phyla Spirochaetes (one genera) and Verrucomicrobia (one genera).
These results were also reported in ruminant research on animals such as Yaks and Tibetan sheep [21,22], indicating that although TWAs and NPDDs belonged to the genus Equus and were different from ruminants, their microorganism flora communities were actually quite similar. In addition, an average of 54.43% genera presented as unclassified, which suggested that this group of organisms may be a part of the core bacterial population, warranting further investigation in the field on Equus.
Recent studies have shown that changes in gut microbiota composition are associated with host disease [5,23,24]. Bacteroidetes, a basal microbe, is reported to be the most abundant phylum in healthy people [25] and captive TWAs [12]; a decrease in its relative abundance is associated with chronic diarrhea in humans [26] and could be a predictor of an animal's health. In our study, it accounted for a large proportion in TWAs (44.12% of sequences), indicating that TWAs may have a superior capacity to stave off intestinal diseases. The phylum Spirochaetes, which contained pathogens that cause intestinal disease [27], were richer in NPDDs than in TWAs, which indicated that NPDDs would more likely suffer from intestinal diseases. At the genus level, Akkermansia accounts for the highest proportion in Verrucomicrobia. Researchers have found that Akkermansia contains a single species, namely Akkermansia muciniphila, which is an important probiotic that contributes healthy mucus-related microbiota composition [28,29]. In addition, Akkermansia muciniphila can be used as a probiotic to prevent obesity and type 2 diabetes [24,30]. In our research, TWAs had the higher distribution proportion of Akkermansia muciniphila, which may be important to the maintenance of intestinal health; however, further study it is necessary to verify this assertion.
The metabolites of gut microbiota play a pivotal role in maintaining physiologic and metabolic homeostasis in their hosts [31]. In the present study, the metagenomic function prediction showed that the two predominant intestinal microbiota in TWAs and NPDDs were related to ABC transporters and the two-component system (Table S4). Previous research [32,33] has confirmed that ABC transporters constitute one of the largest known protein families and are widespread in bacteria, archaea, and eukaryotes, and the twocomponent system is a signal transduction system that senses developmental and environmental stimuli [34]. Therefore, it is reasonable that these KOs were found to be in high abundance in the intestinal microbiota of TWAs and NPDDs. Moreover, statistical comparisons showed that 40 different metabolic pathways exhibited significant differences (Table S2). Among them, 29 pathways had richer concentrations in TWAs than NPDDs, mainly amino acid metabolism, carbohydrate metabolism, and energy metabolism.
A large number of various microorganisms inhabit the intestinal tracts of monogastric animals that are symbiotic with the hosts. Gut bacteria produce a range of metabolites that provide nitrogen sources for themselves and amino acids for the hosts through the fermentation of protein and nitrogen compounds [35]. In our study, the high concentration of amino acid metabolic pathways in TWAs implied that in the context of high altitude and forage shortage, the gut microorganisms and hosts coevolved and exploited high bioavailability to efficiently utilize nitrogen compounds in forage, ultimately with microorganisms providing themselves a nitrogen source and in return producing various amino acids to maintain normal life activities for the hosts. Research has revealed that carbohydrates in forage can be fermented by intestinal microorganisms and the metabolites include certain short chain fatty acids (SCFAs) which can stimulate the intestinal mucosa, improve intestinal immunity, and promote energy absorption in the hosts [36,37]. In our study, the concentration of carbohydrate metabolic pathways in TWAs had an advantage over NPDDs, which meant that microorganisms in the gut of TWAs could more efficiently promote cellulose and hemicellulose fermentation. As for energy metabolism, we found that the citrate cycle in TWAs was significantly more active than in NPDDs, indicating that the intestinal bacteria of TWAs could provide more energy for their own metabolic activities and were therefore beneficial for microorganism fermentation of forage in the intestine.
There are complex interactions in gut microbiota, such as mutualism, competition, prediction and commensalism [38]. Understanding the interactions and their responses to environmental is an important goal in ecology. However, defining a network structure in microbial community is especially challenging due to their vast diversity and as-yet uncultivated status. Thus, we used the network analysis to investigate the genus interactions according to their mutual competition (negative correlation) and mutual benefit (positive correlation). A positive correlation may relate to mutualism, commensalism or parasitism, while a negative correlation may result to competition, predation, etc [38]. In our network analysis, the TWAs were predominated by a higher positive correlations, while NPDDs were predominated by negative correlations. We speculated that plant-based diet in the gut of TWAs was more conductive to fermentation through the more corporation interactions than that of NPDD. In addition, we found that in TWAs, the total links, avgD and D were relatively higher than that of NPDDs, indicating that TWAs harbored a more complex ecological network, mostly also more metabolite.
More complex interactions could build a more stable metabolic network environment for microbial population. Further study is needed to confirm this hypothesis.
A lower avgCC revealed that bacterial communities are largely composed of isolated nodes or loosely connected node groups, reflecting lower functional redundancy [39]. In our study, the avgCC in TWAs was relatively higher than NPDDs, indicating that the latter may harbor lower function redundancy than the former. Previous study had revealed that centralized nodes represent keystone species, which may influence the "information" flow among microbes [40]. In addition, these nodes play an important role in connecting many other microbes and are generally considered as important control points in the networks A low level of CB, as a result of low frequency of centralized nodes [40].Thus, a higher CB in TWAs, indicated that the microorganism had different status in networks, implying that the lose of one or several keystone species may affect the stability of the whole bacterial community. Our result also showed that NPDDs harbored more key stone genus than TWA.
Here, a total of three module hubs and six connectors were observed in TWAs (four connectors) and NPDDs (three module hubs and two connectors, respectively). However, none of them were the top 10 most abundant genus. The most abundant module microorganism OTU67 was at the 53th position, which indicated that the key stone genus in network not necessary the dominant microbes in the gut microbial community.
However, the topological features are related to sequencing depth and different data processing. In our study, we only remained the OTUs that were present more than half of all samples, and this data processing may influence the topological features in bacteria network.

Conclusion
This study described the compositions, diversity and function of bacterial commubities in Tibetan wild assed and domestic donkeys. Based on 16s rRNA sequencing data, we found that the gut microbiota of TWAs was superior to that of NPDDs in bacterial community composition, function, and potentially high resistance to disease risk under similar forage nutrition intake. As the only wild species of perissodactyls on the Qinghai-Tibetan Plateau, the protection of Tibetan wild asses is of great significance to maintain the balance of the ecosystem and the harmony between human and animals. From the perspective of molecular biology, we studied the difference of gut microbiota between Tibetan wild asses and domestic donkeys. These results may help us understand the assembly of bacteria community in the gut of wild and domestic animals, and provide a theoretical basis for Tibetan wild asses' protection.

Experimental design and sample collection
In this study, all animal care procedures were consistent with the guidelines of Sanjiangyuan National Park Administration and Northwest Institute of Plateau Biology, CAS (NWIPB20160302).
For the feces collection of natural pasture domestic donkeys (NPDDs), the study was On the morning of Aug 16 th , each donkey was attached to the grassland with a rope about 10 meters long, and the feeding range was kept about 15 meters apart. We numbered the 15 laboratory animals from 1 to 15, and observed them in the corner where they were not easy to find. Once the fresh feces were discharged, according to the number, we collected them immediately and stored in sterile self-sealing bags with 4 degree portable refrigerator. In this way, all fresh feces that each donkey natural defecation were collected during the day. Finally, mixed the faces sufficiently, transformed into sterilized freezing tubes (approximately 2 milliliters each donkey) and stored in liquid nitrogen for later DNA extraction.
For the feces collection of Tibetan wild asses (TWAs), the study was conducted in alpine meadow (101°.06′E, 38°.37′N) at Qumarleb County (Yushu Prefecture, Qinghai Province, China) with an altitude of 4300 m in August 2017. This is a habitual area for Tibetan wild asses. On the morning of Aug 3 th , we found a herd of 35 TWAs, through binoculars observation of the wild animals' behavior, we selected 15 adult, females and healthy TWAs as the laboratory animals. We used binoculars to observe the defecation of each wild ass, and 2 milliliters sterilized freezing tubes were used to quickly collect the mixed fresh feces and stored in liquid nitrogen. During the day, we walked about 20,000 steps and collected the feces of 15 experimental wild asses. During sampling, all laboratory wild animals had free access to water, no scared, no driven and antibiotic injection.
The forage samples were collected from quadrats (50 cm × 50 cm) of grass on which the animals grazed. Twenty quadrats greater than 10 m apart were placed randomly to investigate the dominant species and collect the ground herbage. A total of 20 forage samples were collected from these two regions (10 per region). Forage samples were dried in a 65℃ oven for 24 h, then ground through a 1-mm sieve and stored in a vacuum dryer for nutritional analysis.

Determination of plant nutrient composition
Total N was measured using the Kjeldahl method; the crude protein content was calculated as 6.25 × N (Method No. 984.13); the ether extract (EE) was measured using the Soxhlet system (Method No. 954.02); the acid detergent fiber (ADF) and neutral detergent fiber (NDF) were analyzed using the method described by Soest et al. (1991) [41]; the non-fibrous carbohydrate (NSC) was calculated as follows: NSC % = 100-Ash%-EE%-NDF%-CP(%). The nutrition of the herbage and dominant species are shown in Table   1.

Sequencing and processing
To get clean paired-end reads, the raw reads were filtered to remove reads containing greater than 10% unknown nucleotides (Ns) and fewer than 80% of base calls with quality scores less than 20. FLASH version 1.2.11 was used to merge paired-end reads as raw tags with a minimum overlap of 10 bp and mismatch error rates of 2% [42]. The noisy raw tag sequences were filtered by QIIME to obtain high-quality, clean tags [43]. The clean tags then underwent chimera detection using the UCHIME algorithm until allchimeras were removed and effective tags were obtained [44]. All of the sequences were clustered into operational taxonomic units (OTUs) of "≥ 97%" similarity using the UPARSE version

Statistical analysis
Taxonomic composition for each cluster was evaluated at the phylum and genus level.
Data related to bacterial community were statistically analyzed using SPSS version 17.0 (SPSS, Inc., Chicago, IL, USA). The significance of herbage nutritional composition and bacterial taxa was determined using the independent-samples T test. The bacterial community function were predicted using Tax4Fun software [48]. The significant KEGG metabolic pathways were screened using STAMP software (http://kiwi.cs.dal.ca/Software/STAMP). The OTU rarefaction curves were calculated and plotted in QIIME. All alpha diversity was calculated in QIIME and graphed by Origin version 8.0 (OriginLab®, Northampton, MA, USA)). The principal coordinates analysis (PCoA) with unweighted UniFrac distance was plotted in R version 3.5.0. The unweighted pair-group method with arithmetic mean (UPGMA) classification tree was plotted by mothur version 1.41.1 [49]. Results are reported as means ± SD. Effects were considered significant at P < 0.05.

NPDDs
In order to understand the interaction of gut microbiota between TWA and NPDD, phylogenetic molecular ecological networks (pMENs) (http://ieg4.rccc.ou.edu/mena) were constructed based on Random Matrix Theory (RMT)-based methods. The protocols of network construction were described previously [50]. Briefly, only the OTUs that were presented in more than half of all samples were selected to in each TWAs and NPDDs gut microbiota genus. In order to compare the topological characteristics of bacteria network between TWAs and NPDDs, the pMENs were calculated with the same threshold (0.76). For each network, in above website, the window "global network properties" was used to calculate total nodes, total links, average degree (avgK), average clustering coefficient (avgCC), centralization of betweenness (CB) and density (D). The window "module separation and modularity calculation" was then used to calculate the value of withinmodule connectivity (Zi) and among modularity connectivity (Pi). Finally, the visualization of these community co-occurrence networks was made with  Values in the same row with different superscripts are significantly different (P < 0.05).    The principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) of metabolite pathways between TWAs and NPDDs.

Figure 6
Comparision of gut bacterial networks between TWAs and NPDDs. Highly positive correlations is indicated by red color and negative correlations by gray color.TWA=Tibetan wild asses, NPDD=natural pasture domestic donkeys.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.