Comparative study of gut microbiota in Tibetan wild asses (Equus kiang) and domestic donkeys (Equus asinus) on the Qinghai-Tibet plateau

Sample collection

A total of 30 fecal samples were obtained from two different sites. The TWAs group included 15 fecal samples from a swamp meadow (101°.06′E, 38°.37′N) in Qumarleb County (Yushu Prefecture, Qinghai Province, China) at an altitude of 4,300 m in August 2017; the NPDD group included 15 fecal samples from an alpine meadow (94°48′E, 34°.99′N) in Menyuan County (Haibei Prefecture, Qinghai Province, China) at an altitude of 3,200 m in August 2017. All of the laboratory animals were healthy and none of the donkeys had been administered antibiotics within the past three months. The fecal samples (approximately 2 ml each) were collected and transferred into separate sterilized tubes and stored immediately in liquid nitrogen for DNA extraction. 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 °C 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 fibre (NDF) were analyzed using the method described by Soest, Robertson & Lewis (1991); the non-fibrous carbohydrate (NSC) was calculated as follows: $\mathrm{NSC}\left(\text{\%}\right)=100-\mathrm{Ash}\left(\text{\%}\right)-\mathrm{EE}\left(\text{\%}\right)-\mathrm{NDF}\left(\text{\%}\right)-\mathrm{CP}\left(\text{\%}\right)$. The nutrition of the herbage and dominant species are shown in Table 1.

Acid-insoluble ash to determine dry matter digestibility

The dry matter digestibility was measured by using a modified method of 4N-HCL acid-insoluble ash (AIA) (Van Keulen & Young, 1997). Briefly, 5 gram crushed forage and 5 gram dried excreta are boiled in 50ml 4N for thirty minutes; the slurry is then filtered by quantitative filtered paper (NEWSTAR®, Hangzhou, China, Φ 12.5 cm) and the residue washed free of acid using hot distilled water untilled the elute is neutral; the residue and the filter paper were ashed directly in crucible at 600 °C for 12 h in electric muffle furnace (Beijing zhongxing weiye instrument co. LTD, KSW-6-12); After the residue and filter paper complete combustion, the crucible and ash were placed in a dryer to cool to room temperature and then weighed on an analytical balance. The equation used to calculate dry matter digestibility was as follows (Kavanagh et al., 2001): Dry matter digestibility = 100 − 100 × AIA in herbage∕AIA in faeces.

DNA extraction and purification

Microbial DNA was extracted from stool samples using the E.Z.N.A® Stool DNA Kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer’s protocols. DNA quality was assessed via 2% agarose gel electrophoresis and metagenomic DNA concentrations were determined with a NanoDrop™ 2000 (Thermo Fisher Scientific, Waltham, MA, USA). The 16S rDNA V3-V4 regions of the ribosomal RNA gene were then amplified by PCR (95 °C for 2 min, followed by 27 cycles at 98 °C for 10 s, 62 °C for 30 s, and 68 °C for 30 s, with a final extension at 68 °C for 10 min) using the following primers: 341F, CCTACGGGNGGCWGCAG; 806R, GGACTACHVGGGTATCTAAT. The 5′-end of 314F primer includes 8-bp unique barcodes, which were used to split each sample. PCR reactions were performed in triplicate with a 50 µl mixture containing 5 µl of 10 × KOD buffer, 5 µl of 2.5 mM dNTPs, 1.5 µl of each primer (5 µm), 1 µl of KOD polymerase, and 100 ng of template DNA. The PCR products were detected by electrophoresis on a 2% agarose gel and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer’s instructions and quantified using QuantiFluor®-ST (Promega, Madison, WI, USA). The purified amplicons were then pooled equimolarly for paired-end sequencing on an Illumina (San Diego, CA, USA) HiSeq 2500 platform (Guangzhou Gene Denovo Biotechnology Co Ltd Guangzhou, China) using standard protocols.

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% (Magoc & Salzberg, 2011). The noisy raw tag sequences were filtered by QIIME (version 1.9.1) to obtain high-quality, clean tags (Caporaso et al., 2010). The filter commands mainly including tags interception and length filtering. That is, the raw tags were intercepted from the continuous low-quality value (the default quality threshold is <=3) to the first low-quality base site of the set length (the default length is 3); filtering out tags data that had a continuous high-quality base length less than 75% of the length of the tags. The clean tags then underwent chimera detection using the UCHIME algorithm until all chimeras were removed and effective tags were obtained (Edgar et al., 2011). All of the sequences were clustered into operational taxonomic units (OTUs) of “ ≥ 97%” similarity using the UPARSE version 7.0.1001 pipeline (Edgar, 2013) and the representative sequences were classified into organisms by applying a naive Bayesian model using the RDP classifier version 2.2 (Wang et al., 2007) based on the SILVA database (version v132 ) (Pruesse et al., 2007).

Statistical analysis

The 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 was predicted using Tax4Fun software (version 0.3.1) (Asshauer et al., 2015). The significant KEGG metabolic pathways were screened using STAMP software (version 2.1.3), the statistical method was two groups Welch’s t-test, the type was two-sided and the confidence interval method was Welch’s inverted (0.95). 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 at OTU level was plotted in R version 3.5.0 and metric distance were calculated by Wilcoxon rank sum test. The unweighted pair-group method with arithmetic means (UPGMA) classification tree was plotted by mothur version 1.41.1 (Schloss et al., 2009). Results are reported as means ± standard error (SE). The effects were considered significant at P < 0.05.

Network analysis and identification of putative keystone taxa in TWAs and NPDDs

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 (Li et al., 2016). 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 the above website, the window “global network properties” was used to calculate total nodes, total links, positive links, negative links, average degree (avgK), average clustering coefficient (avgCC), centralization of betweenness (CB) and density (D). It was worth noting that a positive correlation might relate to mutualism, commensalism or parasitism, while a negative correlation might result to competition, predation, etc. (Faust & Raes, 2012). The window “module separation and modularity calculation” was then used to calculate the value of within-module connectivity (Zi) and among modularity connectivity (Pi). Finally, the visualization of these community co-occurrence networks was made with Cytoscape 3.3.0 (Shannon et al., 2003).

A module typically contains many nodes, which are tightly connected to each other in one group, while having only few connections outside the group. The values of Zi and Pi are indicators of the connectivity of each node, and thus are often used to determine the topological role of these nodes (Deng et al., 2012). According to the valises of Zi and Pi, these nodes can be classified into four categories (Zhou et al., 2011), including peripherals (In the modules, the OTUs have few outside connections, Pi < 0.62 and Zi < 2.5), connectors (OTUs that connect modules, Pi > 0.62), network hubs (OTUs that highly connected within entire network, Pi > 0.62 and Zi > 2.5) and module hubs (OTUs that highly connected within modules, Zi > 2.5). Connectors, network hubs and module hubs were considered as the putative keystone taxa in a microbial community (Li et al., 2017). The keystone of gut microbiota in TWAs and NPDDs was identified based on the above values.

Analysis of food nutrition composition and dry matter digestibility between TWAs and NPDDs

As shown in Table 1, the mainly plant-based diet for TWAs were Kobresia tibetica and Kobresia pygmaea, while for NPDDs were Kobresia humilis and Kobresia capillifolia. There was no significant difference in nutrient composition except ash content (P < 0.05). In addition, the dry matter digestibility was measured by acid-insoluble ash method. As shown in Table 2, AIA in the faces of NPDDs were significantly higher than TWAs group (P < 0.05), whereas, dry matter digestibility in NPDDs were significantly lower than TWAs group (P < 0.05).

Sequencing and classification

Illumina sequencing yielded 2,540,672 raw reads of 16S gene sequences. After quality filtering, a total of 2,395,867 effective tags were obtained (Table S1). As shown in Fig. S1, 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 Tibetan wild ass and domestic donkey

To determine alpha diversity, we calculated the Chao1, ACE, observed species estimates and the Shannon diversity index. As shown in Fig. 1 and Fig. S2, no significant differences were found in these two groups for the Shannon index, Chao1, ACE estimates, the number of observed species and Faith’s phylogenetic diversity (PD). The OTU level rarefaction curves of diversity estimators reached a plateau at the minimum sequence number of 50,000 tags (Fig. S3). With regard to beta diversity, PCoA analysis showed that the two groups of bacterial communities differed significantly (Fig. S4A) (P = 3.02E-05). To assess the overall similarity, UPGAMA with classification tree analysis was conducted as shown in Fig. S4B. 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 and Table S2. 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, Spirochaetaes, Proteobacteria, and Planctomycetes between the two groups (P > 0.05).

At the genus level, a total of 15 genera with proportions above 1% were detected, as shown in Fig. 3A. The Rickenllaceae_RC9_gut group, Fibrobacter, and Treponema_2 accounted for the largest proportions in all samples. Differentiation analysis of the 15 genera is presented in Fig. 3B and Table S3. The proportion of Ruminococcaceae_NK4A214, Phascolarctobacterium, Coprostanoligenes_group, Lachnospiraceae_XPB1014_group, and Akkermansia in TWAs was significantly higher than in NPDDs, whereas the proportion of the Lachnospiraceae_AC2044_group in NPDDs was significantly higher than in TWAs (P < 0.05). No remarkable difference was detected in the relative abundance of the Rickenellaceae_RC9_gut_group, Fibrobacter, Treponema, Ruminococcaceae_UCG-010, Christensenellaceae_R-7_group, Prevotellaceae_UCG-001, Anawrovorax, Prevotellaceae_UCG-003, and Prevotellaceae_UCG-004 between the two groups (P > 0.05).

To study the main differences in gut microbiota between groups TWA and NPDD, linear discriminant analysis (LDA) combined effect size measurements (LEfSe) were used to identify the discriminative features of bacteria in the TWAs and NPDDs (Fig. S5). The taxonomic units of Campylobacter, Campylobacteraceae, Campylobacterales, Epsilonproteobacteria, and Campylobacter_hyointestinalis_subsp_hyointestinalis, which belong to the phylum Proteobacteria, were significantly enriched and revealed different biomarkers in TWAs, whereas no specific taxonomic units were enriched in NPDD.

Predicted metabolic pathways and functions of the bacterial microbiota between TWAs and NPDDs

A total of six types of biological functional pathways in KEGG level 1 were detected (Fig. 4A). The metabolism pathways had the highest relative abundance, with more than 60% of the total reads in each group. At KEGG level 2, a total of 37 metabolic pathways were detected in the gut samples (Table S4). 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 S5, 40 metabolite pathways were significantly different (VIP > 1 and P < 0.05), of which ten metabolite 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 metabolites were 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 the 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

As shown in Fig. 6, when using the same threshold (0.76), a total of 207 and 104 links were identified for TWA and NPDD. Notably, for TWAs, the percentage of positive correlations (149 out of 207) was dominant, while negative correlations (55 out of 104) were the predominant in NPDDs. Interestingly, we found the values of total links, avgK and D were higher in TWAs than those in NPDDs, indicating a relatively simpler network in NPDDs. The lower avgCC in NPDD indicated that the bacterial network was mainly loose node groups. Moreover, TWA had a higher level of CB than NPDD, implying that TWA harbored a higher frequency of centralized.

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