The Composition and Predictive Function of the Fecal Microbiota Differ Between Young and Adult Donkeys

The community of microorganisms inhabiting the gastrointestinal tract of monogastric herbivores played critical roles in the absorption of nutrients and keeping the host healthy. However, its establishment at different age groups has not been quantitatively and functionally examined. The knowledge of microbial colonization and its function in the intestinal tract of different-age donkeys is still limited. By applying the V3–V4 region of the bacterial 16S rRNA gene and functional prediction on fecal samples from different-age donkeys, we characterized the gut microbiota during the different age groups. In contrast to the adult donkeys, the gut microbiota diversity and richness of the young donkeys showed significantly less resemblance. The microbial data showed that diversity and richness increased with age, but a highly individual variation of microbial composition was observed at month 1. Principal coordinate analysis (PCoA) revealed a significant difference across five time points in the feces. The abundance of Bacteroides, Lactobacillus, and Odoribacter tended to decrease, while the proportion of Streptococcus was significantly increased with age. For functional prediction, the relative abundance of pathways had a significant difference in the feces across different age groups, for example, Terpenoids and Polyketides and Folding, Sorting, and Degradation (P < 0.05 or P < 0.01). The analysis of beta diversity (PCoA and LEfSe) and microbial functions predicted with PICRUSt (NSTIs) clearly divided the donkeys into foals (≤3 months old) and adults (≥7 months old). Microbial community composition and structure had distinctive features at each age group, in accordance with functional stability of the microbiota. Our findings established a framework for understanding the composition and function of the fecal microbiota to differ between young and adult donkeys.

Introduction

Monogastric herbivores are an important component of the agricultural sector, due to the meat, milk, and skin they produce for human use; however, these animals must consume increasing feed resources to meet the demands of the growing human population (Eisler et al., 2014). Consequently, improvement of feed efficiency is critical to the development of a sustainable monogastric herbivore. In particular, the major end products of the gut microbiota catabolism heavily depend on dietary fiber, i.e., the short-chain fatty acids acetate, propionate, and butyrate; these are a key energy source for the monogastric herbivores (Brokner et al., 2016). Indeed, the gut microbial systems, which aid in the digestion of otherwise indigestible nutrients, are the most complex biological system in animals and have been considered as an extra digestive organ (Eisler et al., 2014). A growing body of evidence demonstrates that the relatively stable function of the gut microbiota supplies the host with increased adaptability to the environment (Gao et al., 2020) and plays fundamental roles in intestinal physiological development and nutrient metabolism (Bäckhed et al., 2015; Wu et al., 2016). Therefore, understanding the gut microbiota of donkeys at different age in-depth not only can increase animal health, feed utilization efficiency, and production, but also may provide guidelines to early weaning technique.

Donkeys and horses (genus Equus) are typical monogastric herbivores. Almost all herbivores use crude fiber as a fermentation substrate in the intestine, which is the main activity location of microbiota (Bergman, 1990). Many diseases in the horse are related to the destruction of the intestinal microbiota, the most obvious being colitis (Harlow et al., 2013), laminitis (Moreau et al., 2014), and transient diarrhea (Kuhl et al., 2011) in young foals. At present, the gut microbiota of the newborn equine is already complex; for example, Satokari et al. (2008) isolated Streptococcus and Enterococcus from the umbilical cord and meconium of newborn foals. It has been shown that the delivery mode and breast-feeding strongly affect the development and function of the intestinal microbiota in humans (Bäckhed et al., 2015). The initial microbial colonization alters the intestinal environment to benefit the growth of specific anaerobes, which plays a dominant role in later growth and development (Favier et al., 2002). The gut microbiota thus clearly plays an important role for equine health and forage utilization. At present, most studies of equine intestinal microbiota focus on feces (Liu et al., 2014), different regions of the digestive tract (Schoster et al., 2017), and the characteristics of the microbiota in healthy and diarrheic animals (Costa et al., 2012). From an ecological point of view, early colonization of the gut microbiota represents the de novo assembly of a microbial community (Costello et al., 2012) and is influenced by animal species, age, and diet (Zoetendal et al., 2001; Torok et al., 2008; Carmody et al., 2015).

Donkeys are an important species as its skin could be used to produce Ejiao, a famous Chinese medicine, and are also considered as highly valued livestock due to their high-quality meat and skin. Now, the Chinese livestock industry includes approximately 2,500,000 donkeys. In previous studies, we have examined the different digestive tract microbiota of donkeys (Liu et al., 2019). However, limited information is known about the development of the bacterial community composition and function in the feces of donkeys at different ages. Therefore, in this study, we aimed to determine how microbial diversity and function changes in donkeys’ feces with age. The results of this study may shed new light on the microbiota–functional interactions at different ages.

Materials and Methods

Animals, Management, and Sample Collection

All procedures involving animals were approved and authorized by the Animal Welfare Committee of Liaocheng University. This study included 25 healthy DeZhou donkeys: they were divided into five groups (F1, F3, F7, F12, and F24) according to their age, and five heads in each group. Donkeys were kept outdoors at the Dong-E-E-Jiao Co., Ltd., national black donkey-breeding center (Donge county, Shandong Province) (June 25, 2019, temperature 23–34°C) and suckled until weaning (month 7). All donkeys had no gastrointestinal diseases or any antimicrobial exposure during the previous 3 months. They were fed a standard concentrate diet comprising 1.3% body weight twice daily (at 8:00 a.m. and 4:00 p.m.), and soybean straw (ratio of 60:40) had a free feeding (Burden and Bell, 2019). Meanwhile, they had free access to fresh water. Unweaned donkeys and their mother donkeys havea free access to soybean straw, while only unweaned donkeys could intake the concentration in their stables. Foals were weaned when they were at month 7 by physical separation from their mother.

The farm technician would collect fecal samples from donkeys after feeding 1 h at months 1, 3, 7, 12, and 24 after birth, picking feces out of the rectum with sterile gloves. Samples were immediately sealed in 50-ml aseptic cryopreservation tubes and put in liquid nitrogen and then stored at −80°C until DNA extraction.

Bacterial DNA Extraction, Purification, and High-Throughput Sequencing

We extracted total genomic DNA from each fecal sample animal using a QIAam stool Mini Kit (Qiagen, Valencia, CA, United States), according to the manufacturer’s instructions. The purity and integrity of the DNA were checked using gel electrophoresis.

The bacterial primer pair 341F (5′-CCTAYGGGRBGC ASCAG-3′) and 806R (5′-GGACTACNNGGG TATCTAAT-3′) (Sonnenburg et al., 2010) were used to amplify the V3–V4 regions of the bacterial 16S rRNA gene (Liu et al., 2019) using Phusion High-Fidelity PCR Master Mix (NEB, New England Biolabs). Each primer contained the appropriate Illumina adapter sequence and 8-bp identifier indices. The resulting amplicons were purified using QIA quick Gel Extraction Kits (Qiagen, Valencia, CA, United States), and pooled in equimolar concentrations. Then, the purified amplicons were quantified using a fluorescence spectrophotometer (Thermo Scientific, MA, United States). Constructed libraries were subjected to Qubit quantitation and library testing, and then sequenced on an Ion S5TMXL system (Thermo Fisher Scientific, Franklin, MA, United States).

Bioinformatics Analysis

Raw sequencing data typically contains a certain proportion of low-quality sequences. To ensure the accuracy and reliability of follow-up analyses, the original data (raw reads) must be quality controlled and filtered to obtain effective data (clean reads) (Haas et al., 2011). We removed low-quality reads using Cutadapt (v1.9.1) and barcode and primer sequences using Chromas. Sequences were spliced together using DNASTAR and removed chimeric sequences using Bellerophon (Huber et al., 2004).

The cleaned sequences were analyzed using QIIME (v1.7.0) (Caporaso et al., 2010) and clustered into operational taxonomic units (OTUs) based on 97% identity using UPARSE (Edgar, 2013). The representative sequences of each OTUs were aligned to the Greengenes database (DeSantis et al., 2006) using PyNast (Gregory et al., 2009). After that, the taxonomic information was annotated against the SILVA database using the RDP Classifier with a 0.80 confidence threshold (Wang et al., 2007; Christian et al., 2012). OTUs abundance information was normalized using a standard of sequence number corresponding to the sample with the least sequences.

To identify shared OTUs among groups, we constructed a Venn diagram using the Draw Venn Diagram online tool (Heberle et al., 2015). Alpha diversity indexes (i.e., Good’s coverage, Chao1, Ace, observed species, and Shannon) using QIIME (v 1.7.0) (Gregory et al., 2009) were used to compare the bacterial richness and diversity among groups. We identified statistically significant differences in alpha diversity, and microbial composition among donkeys of different ages was evaluated by one-way AVOVA. Beta diversity (i.e., differences in species complexity among samples) was determined based on principal coordinate analysis (PCoA) using weighted_unifrac and linear discriminant analysis (LDA) effect size (LEfSe); these analyses were performed using Microbiome Analyst, a comprehensive web server for comparative metagenomics (Achal et al., 2017). Community membership and structure were compared and visualized between pairs of samples using the PCoA. LEfSe was performed to identify the bacterial taxa that differed between pairs of groups (Segata et al., 2011). Taxa with an LDA score > 4 were considered important biomarkers of each group. We considered P value < 0.05 to indicate significant differences among groups.

Predicted Function of the Gut Microbiome

To reveal the potential metabolic capabilities of the bacterial communities, we used PICRUSt (phylogenetic investigation of communities by reconstruction of unobserved states) to predict the function of gut microbiota based on the OTU table (Langille et al., 2013). Before PICRUSt analysis, we used the closed-reference OTU-picking protocol in QIIME against Greengenes to construct PICRUSt-compatible OTUs. Then, PICRUSt was used to predict functional pathways based on Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations (Breiman, 2001; Segata et al., 2011) at levels 1 and 2 (Langille et al., 2013). We used the nearest sequenced taxon index (NSTI) to evaluate the accuracy of the predictions of the metagenomes. This algorithm used a phylogenetic tree of 16S rRNA gene sequences to associate OTUs with gene content. The relationships between functional capacities and predicted relative gene abundances were analyzed using PCoA and heatmaps. LSD-t was used to compare the relative abundance changes between different groups.

Results

Sequence Data and Alpha Diversity

We obtained 1,655,064 16S rRNA gene sequences across all fecal samples. The sequences were clustered into 3,235 OTUs at 97% sequence identity: 26 phyla, 33 classes, 62 orders, 105 families, and 209 genera. Core OTUs are common OTUs in each sample. We identified 481 core OTUs in all samples and 448 OTUs unique to one age group: 80 OTUs unique to group F1, 79 OTUs unique to group F3, 106 OTUs unique to group F7, 101 OTUs unique to group F12, and 82 OTUs unique to group F24 (Figure 1A). The rarefaction curves across all samples tended to plateau at a sequencing depth of 45,802, indicating that this sequencing depth captured most of the microorganisms presented in the samples (Figure 1B). However, not only fewer OTUs consisted in the younger donkeys (groups F1 and F3) (P < 0.01), but lower sequencing depths were also required to cover the fecal microbiota at these ages (Figure 1B).

FIGURE 1

Figure 1. Numbers of bacterial operational taxon units (OTUs) in the feces of donkeys at different ages. Venn diagram, where areas of overlap indicate the numbers of OTUs shared among the overlapping groups (A). Rarefaction curve showing the number of OTUs recovered at different sequencing depths (B). F1: at month 1; F3: at month 3; F7: at month 7; F12: at month 12; F24: at month 24.

There were no significant differences in the numbers of effective sequences obtained across the 25 fecal samples (P > 0.05; Table 1). At a sequencing depth of 45,802, the Good’s coverage indexes ranged from 0.994 to 0.998, indicating that more than 99% of all bacterial taxa were captured across all samples. The species richness indexes (i.e., observed species, Chao1, and ACE) and the species diversity index (Shannon) for the older groups (F7, F12, and F24) were significantly higher than those in the F1 and F3 groups (P < 0.01). There were no significant differences in any indexes among older groups (P > 0.05). In addition, the Shannon index of F1 group was the lowest (4.43), indicating that the fecal microbial diversity of the F1 group was low and the microbiota was relatively simple.

TABLE 1

Table 1. Sequencing statistics for the fecal microbiota of donkeys at different ages, including numbers of OTUs (97% similarity) and diversity indexes.

Taxonomic Composition of Feces at Different Ages

We identified the composition of the top 10 bacterial phyla in the fecal samples (Figure 2A). The phyla Firmicutes and Bacteroidetes were the most abundant, together accounting for about 90% of the OTUs identified across all groups. Firmicutes and Bacteroidetes comprised 54.79 and 39.38%, respectively, of the bacteria in group F1; 50.54 and 40.13%, respectively, of the bacteria in group F3; 46.45 and 45.64%, respectively, of the bacteria in group F7; 55.25 and 35.78%, respectively, of the bacteria in group F12; and 54.88 and 37.86%, respectively, of the bacteria in group F24 (Figure 2A). Although genus with relative abundances less than 1% were most common in all age groups, group F1 had fewer of these genera (35.50%), and much higher relative abundances of Bacteroides (17.34%), Odoribacter (17.20%), Lactobacillus (12.31%), and unidentified Christensenellaceae (8.59%) (Figure 2B). In groups F3–F24, genera with relative abundances less than 1% accounted for 80.79, 86.32, 79.51, and 76.90% of all bacteria, respectively, with the remaining bacteria primarily Streptococcus, unidentified Ruminococcaceae, and unidentified Clostridiales. Interestingly, group F3 retained some Bacteroides (3.45%) and Lactobacillus (2.19%), which were not present in the older samples. Also, the relative abundance of Streptococcus increased with age. Then, we compared the relative abundance of these taxa associated with the different ages (Supplementary Figure 1). The proportion of Bacteroides, Lactobacillus, and Odoribacter reduced linearly from groups F1 to 24 (P < 0.05 or P < 0.01, Supplementary Figure 1B). The abundance of Streptococcus was significantly higher at groups F12 and F24 than at groups F1 and F3 (P < 0.01), and there was no significant difference with group F7 (P > 0.05).

FIGURE 2

Figure 2. Microbial community composition in the feces of donkeys at different ages. Taxonomic composition at the phylum (A) and genus levels (B).

PCoA revealed that the microbial communities from feces across different ages are varied (Figure 3A). The unweighted unifrac distance analysis showed that the community composition of the feces at groups 1 and 3 was significantly separated from other groups (at groups 7, 12, and 24) (AMOVA < 0.01); PC1 explained 38.46% of the variation among groups.

FIGURE 3

Figure 3. Alterations in fecal microbiota of donkeys at different ages. Principal coordinate analysis (PCoA) showing phylum-level differences (A). Histogram showing taxa with statistically significantly linear discriminant analysis (LDA) values of donkeys at different ages (B). f, family; g, genus; s, species; o, order.

We employed a LEfSe analysis to identify the taxa that most strongly differentiated the age groups. More differential bacteria were found through comparing to other groups. At group 1, there were several microbial populations that were enriched; these included Bacteroidaceae, Marinifilaceae, Odoribacter, and Lactobacillaceae. Only one family was found at groups 3 (Ruminococcaceae) and 7 (Rikenellaceae), respectively. In comparing the different bacteria of fecal samples at groups 12 and 24, there were proportional differences between Prevotellaceae and Spirochetes at group 12, and Streptococcaceae, Lactobacillaceae, and Bacilli at group 24 (Figure 3B).

Functional Predictions of the Rectal Microbiota Using PICRUSt

We predicted the function of the fecal microbiota of donkeys at different ages using PICRUSt, and we used the NSTIs to evaluate the accuracy of these predictions. Metabolism was the most overrepresented level 1 pathway across all groups, while membrane transport, carbohydrate metabolism, amino acid metabolism, and replication and repair were among the most overrepresented level 2 pathways (Figures 4A,B). However, of the level 1 pathways, environmental information processing was significantly more common in group F1 as compared to groups F3 and F7 (P < 0.05) (Figure 4C).

FIGURE 4

Figure 4. Relative enrichment of KEGG Level 1 (A) and Level 2 (B) pathways in fecal microbiota of donkeys at different ages. Corresponding heatmaps for KEGG Level 1 pathways (C) and KEGG Level 2 pathways (D). Differences in relative abundance of each functional gene were examined using one-way ANOVA, and letters (a, b, c, d, e, f, and g) indicate significant differences between groups 1 and 3; between groups 1 and 7; between groups 1 and 12; between groups 1 and 24; between groups 3 and 7; between groups 3 and 12; and between groups 3 and 24, respectively. Lowercase letters correspond to significant differences (P < 0.05); capital letters correspond to extremely significant differences (P < 0.01).

Of the 35 level 2 KEGG pathways, the relative abundance of Membrane Transport, Carbohydrate Metabolism, and Amino Acid Metabolism were at the highest level during different ages. The result shows that metabolism of Terpenoids and Polyketides and Folding, Sorting, and Degradation were significantly gradually enriched from group F1 to group F24 (P < 0.05 or P < 0.01), while Poorly Characterized and Transcription pathways were significantly gradually reduced from group F1 to group F24 (P < 0.05; Figure 4D). Cancers, Nervous system, and Infectious diseases pathways were significantly enriched in group F3 as compared to group F12 (P < 0.05). The 35 level 2 KEGG pathways among F7, F12, and F24 had no significant difference (P > 0.05), which indicated that the prediction function of microorganisms in rectal contents at 7 months is similar to that of adult donkeys, and 7 months of age is a transitional stage. PCoA revealed that the functions of the fecal microbiota from different ages were clustered, with the first two components explaining a total of 67.93% of the variation. The results suggested that the functional KOs at group F1 were distinctly separated from other groups, explaining 41.63% variation (Figure 5).

FIGURE 5

Figure 5. PCoA of KEGG Level 2 pathways in fecal microbiota of donkeys at different ages.

Discussion

In the present study, we described the early colonization process of the donkey gut microbial ecosystems, on the basis of different age. When searching for bacterial OTUs shared among all groups, we were successful in detecting 481 OTUs shared, and each group has its own unique OTUs. The existence of OTUs specifically shared among all groups suggests that each of them contained similar bacterial community, and there was also specificity in terms of microbial DNA components. We also found that microbial diversity in the feces from the youngest donkeys (month 1) was highly variable as compared to those of older donkeys. The species richness and Shannon indexes were higher at 7, 12, and 24 months than at 1 or 3 months, indicating an increase in richness with developmental age. This result was consistent with previous investigations of the feces, ilea, or rumens of foals (Costa et al., 2016), goats (Jiao et al., 2016), deer (Li et al., 2018), and cattle (Uyeno et al., 2010; Georgios et al., 2013) at different ages. Previous studies have indicated that microbial diversity in the gut may be influenced by many factors, including the characteristics of the mother’s milk (in the case of suckling animals) (Mueller et al., 2015), the forage eaten (Harlow et al., 2016), sex (Mshelia et al., 2018), age (Stephens et al., 2016), and the microbiota of the surrounding environment (Chen et al., 2018; Zhou et al., 2019). It has been shown that the bacterial communities in human breast milk are similar to the infant gut microbiota within a few days of birth (Williams et al., 2019), and infants can receive 27.7% of the bacteria from breast milk during the first 30 days of life after receiving 75% of total daily milk intake as breast milk (Collado et al., 2007). However, the components of mother’s milk play a more important role in the early colonization of the intestine with the start of mother’s milk (Quercia et al., 2019), including protein, lactose, milk fat, and other nutrients. Newborn donkeys only have a single stomach and can be physiologically and functionally considered non-ruminants, so the mother’s milk can pass directly to the small intestine. Perhaps this is one of the reasons why microorganisms in the feces of young animals differ from those of adults. The impact of mother’s milk on microbial composition requires further investigation.

The phyla Firmicutes and Bacteroidetes are dominating in the feces over the development of these donkeys. This is consistent with a previous study of the feces of DeZhou donkeys (Liu et al., 2014) and disagrees with the results from the feces of 2-week-old healthy foals (Schoster et al., 2017). This incongruence suggests that the relative abundances of certain bacterial phyla may differ among the feces of animals in the same genus. The Enterococcus and Enterobacteriaceae were found in meconium from birth to the third day of life, indicating that the foal gut microbiota can acquire microorganisms typical of the milk community (Quercia et al., 2019). These results imply that the mother’s milk is an important factor to affect intestine microbial composition. This experiment showed that Bacteroides, Lactobacillus, and Odoribacter are the dominant bacteria in young donkeys (months 1 and 3), and these bacterial genera trend to decrease with age, while Streptococcus is the most abundant bacterial genus in adult donkeys (months 7, 12, and 24). The bacterial composition in feces of different ages was consistent with LEfSe analysis results. The different composition of the genus may be influenced by the mother’s milk while reinforcing the idea that the early microbiota in the feces is related to the milk and the environment. As crude fiber intake increased with age, Streptococcus became the dominant bacterial genus, presumably to break down large volumes of cellulose. However, these hypotheses need to be documented in a future study.

According to our findings, the donkeys’ gut microbial ecosystems describe a specific developmental trajectory from infancy to adulthood, progressively approaching the configuration typical of the adult gut, especially starting from 1 month old to 7 month old. The temporal changes in the gut microbiota of donkeys characterized here indicated that the intestinal flora develops as the donkey ages. In Marinifilaceae and Odoribacter, which were found to dominate at month 1, the bacteria characteristics of the feces were suddenly lost with the start of mother’s milk intake. The high abundance of subsequently uncommon bacterial genera at month 1 might be due to environmental exposure or to the wide range of bacteria to which the foal is exposed through the mother. In infants, the first change of fecal bacterial composition occurs around 5 days of age; it will transform to an adult-like pattern and usually after eating solid food (Bäckhed et al., 2015). In contrast to human, donkeys start to intake small amounts of concentrates already around 10 days of age (Yatsunenko et al., 2012). To our knowledge, the donkeys intake more and more concentrate and roughage with age, so the gut microorganisms of donkeys tend to be more stable with age and more resistant to diseases. The results showed that these taxa do not represent the true colonizers until by month 7, but only transient organisms. The supposition is supported by the quick reduction in the relative abundances of many of the initial microbial genera, as the true colonizers begin to take hold. Our results indicated that month 7 represents a critical period in the intestinal microbiota development, because, by this point, the structure of the intestinal microbiota was similar to that of an adult donkey. From month 7, the donkey gut microbial ecosystems begin to converge to the adult one. This process involves the loss of special microorganisms and the concomitant acquisition of fiber fermenters typical of adult core gut microbiota, such as Prevotellaceae, Spirochetes, Streptococcaceae, and Lactobacillaces. These data indicated that, by month 7, the donkey’s gut microbiota has begun to build up gradually to include the bacteria necessary for the digestion of the roughage found in the adult donkey diet, much like what had been shown in a study of weaning (day 60) and adult (mare) horse (De La Torre et al., 2019). Meanwhile, it is feasible to choose weaning at the age of 7 months in production practice, which not only can ensure the stability of intestinal microorganism but also can reduce the incidence of disease.

PICRUSt functional prediction analysis indicated that, in donkeys of all ages, the membrane transport, carbohydrate metabolism, and amino acid metabolism at level 2 were the most overrepresented. It has been discovered that these pathways exactly were the most active genes in the intestinal tract (Liu et al., 2019). We found that environmental information processing and Poorly Characterized and Transcription pathways were more enriched in the young donkeys at months 1 and 3, while metabolic pathways (i.e., amino acid, energy, enzyme families, and vitamins) were more enriched in the adult donkeys at months 12 and 24. The old donkeys require more energy to survive, so the function of metabolic pathways was enhanced. Also, Cancers and Nervous system pathways were enriched in group F3, and we speculate that the increased intake of forage caused a disturbance in the gut microbiome that stimulated the immunity system of donkeys and drove the inflammation. The differences of microbial function at different ages indicated that the varied pathways enriched in the young donkeys might help to adapt to the post-birth environment and that microbial functions are known to vary throughout the host’s life, due to the interactions between the host and the microbes (Newell and Douglas, 2014). This was consistent with our results; meanwhile, the 35 level 2 KEGG pathways among months 7, 12, and 24 had no significant difference, which suggested that donkeys’ age correlated with differences in the composition, diversity, and function of the fecal microbiota due to the adaptations of the microbiota to their rapidly changing environment. Here, the functional KOs at month 1 was distinctly separated from other groups, indicating that there was a significant difference at month 1 compared with other months and suggesting a shift toward a more adult-like intestinal environment associated with the increased functional capacity for metabolic pathways. Collectively, these results indicate that both the composition and the function of the gut microbiota evolve as the donkeys grow, in response to environment factor and resources feeding into the community. The limitation of the present study is that the metagenomic shotgun sequencing of the feces was not examined, which will provide more accurate and direct evidence to the function of the gut microbiota in future studies.

Conclusion

In summary, we analyzed the fecal microbiota of donkeys at different ages using high-throughput 16S rRNA bacterial gene sequencing and PICRUSt. Our results showed that the microbial diversity in the feces increased as the donkey aged, indicating that age was a key factor in microbial succession. Microbial community composition and functional prediction in the feces of young donkeys (groups F1 and F3) differed significantly from those of older donkeys (groups F7, F12, and F24), implying that the component of mother’s milk influenced the composition of the initial gut microbiota. The results of the five time points revealed that the metabolic pathways changed with age and were more enriched in the adult donkeys at months 12 and 24. Meanwhile, month 7 represented a critical period in the intestinal microbiota development; from month 7, the donkeys’ gut microbial ecosystems began to converge to the adult one. These results increased our understanding of the age-specific variations in the fecal microbiota of donkeys, but additional studies are required to explore the relationships among microbial changes at different ages.

Data Availability Statement

The datasets generated for this study can be found in NCBI SRA accession PRJNA673698, https://www.ncbi.nlm.nih.gov/bioproject/PRJNA673698/.

Ethics Statement

The animal study was reviewed and approved by the Animal Welfare Committee of Liaocheng University.

Author Contributions

JX and GL: laboratory work and writing of the manuscript. SS and DB: experimental development. JY and LL: providing donkeys and sample collection. XW, SS, and YZ: bioinformatic and statistical analyses. GB and MD: study design, interpretation of data, critical review, and final approval of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

The donkey innovation team of Shandong modern agricultural industry technology system (SDAIT-27), the National Natural Science Foundation of China (No. 31860636), and the Agricultural Improved Seed Project of Shandong Province (2017LZGC020) provided funds for experimental animals and results analysis.

Conflict of Interest

JY was employed by the company Dong-E-E-Jiao Co., Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The authors declare that this study received samples from Dong-E-E-Jiao Co., Ltd. Dong-E-E-Jiao Co., Ltd was not involved in the study design, collection, analysis, interpretation of data, the writing of this article or the decision to submit it for publication.

Acknowledgments

We would like to thank Dong-E-E-Jiao Co., Ltd., for providing the donkeys, feeding them, and allowing the experiments to be performed in their farm. We thank their workers for organizing the sample collection.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.596394/full#supplementary-material

Supplementary Figure 1 | Heatmaps of phylum (A) and genus (B) levels abundance. Differences in relative abundance of microbiota at genus level were examined using one-way ANOVA. Letters (a, b, c, d, e, and f) indicate significant differences between groups 1 and 3; between groups 1 and 7; between groups 1 and 12; between groups 1 and 24; between groups 3 and 12; between groups 3 and 24 respectively. Lowercase letters correspond to significant differences (P < 0.05); capital letters correspond to extremely significant differences (P < 0.01).

References

Achal, D., Chong, J., Salam, H., King, I. L., Agellon, L. B., and Xia, J. (2017). Microbiome analyst: a web-based tool for comprehensive statistical, visual and meta-analysis of microbiome data. Nucleic Acids Res. 45, W180–W188. doi: 10.1093/nar/gkx295

PubMed Abstract | CrossRef Full Text | Google Scholar

Bäckhed, F., Roswall, J., Peng, Y., Feng, Q., Jia, H., Kovatcheva-Datchary, P., et al. (2015). Dynamics and stabilization of the human gut microbiome during the first year of life. Cell Host Microb. 17, 690–703. doi: 10.1016/j.chom.2015.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergman, E. (1990). Energy contributions of volatile fatty acids from the gastrointestinal tract in various species. Physiol. Rev. 70, 567–590. doi: 10.1007/BF00370249

PubMed Abstract | CrossRef Full Text | Google Scholar

Breiman, L. (2001). Random forests. Mach. Learn. 45, 5–32.

Google Scholar

Brokner, C., Austbo, D., Naesset, J. A., Blache, D., Bach Knudsen, K. E., and Tauson, A. H. (2016). Metabolic response to dietary fibre composition in horses. Animal 10, 1155–1163. doi: 10.1017/S1751731115003006

PubMed Abstract | CrossRef Full Text | Google Scholar

Burden, F. A., and Bell, N. (2019). Donkey nutrition and malnutrition. Vet. Clin. N. Am. Equine Pract. 35, 469–479. doi: 10.1016/j.cveq.2019.08.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Caporaso, J. G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F. D., Costello, E. K., et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7:335.

Google Scholar

Carmody, R. N., Gerber, G. K., Luevano, J. M., Gatti, D. M., Somes, L., Svenson, K. L., et al. (2015). Diet dominates host genotype in shaping the murine gut microbiota. Cell Host Microb. 17, 72–84. doi: 10.1016/j.chom.2014.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Kaiqian, D., Chao, S., Arunprasanna, V., Xili, L., Yong, L., et al. (2018). Gut bacterial and fungal communities of the domesticated silkworm (Bombyxmori) and wild mulberry-feeding relatives. ISME J. 12, 2252–2262. doi: 10.1038/s41396-018-0174-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Christian, Q., Elmar, P., Pelin, Y., Jan, G., Timmy, S., Pablo, Y., et al. (2012). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596. doi: 10.1093/nar/gks1219

PubMed Abstract | CrossRef Full Text | Google Scholar

Collado, M. C., Derrien, M., Isolauri, E., de Vos, W. M., and Salminen, S. (2007). Intestinal integrity and akkermansiamuciniphila, a mucin-degrading member of the intestinal microbiota present in infants, adults, and the elderly. Appl. Environ. Microbiol. 73, 7767–7770. doi: 10.1128/aem.01477-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. C., Arroyo, L. G., Allen-Vercoe, E., Stämpfli, H. R., Kim, P. T., Sturgeon, A., et al. (2012). Comparison of the fecal microbiota of healthy horses and horses with colitis by high throughput sequencing of the V3-V5 region of the 16S rRNA gene. PLoS One 7:e41484. doi: 10.1371/journal.pone.0041484

PubMed Abstract | CrossRef Full Text | Google Scholar

Costa, M. C., Stampfli, H. R., Allen-Vercoe, E., and Weese, J. S. (2016). Development of the fecal microbiota in foals. Equine Vet. J. 48, 681–688. doi: 10.1111/evj.12532

PubMed Abstract | CrossRef Full Text | Google Scholar

Costello, E. K., Stagaman, K., Dethlefsen, L., Bohannan, B. J. M., and Relman, D. A. (2012). The application of ecological theory toward an understanding of the human microbiome. Science 336, 1255–1262. doi: 10.1126/science.1224203

PubMed Abstract | CrossRef Full Text | Google Scholar

De La Torre, U., Henderson, J. D., Furtado, K. L., Pedroja, M., Elenamarie, O., Mora, A., et al. (2019). Utilizing the fecal microbiota to understand foal gut transitions from birth to weaning. PLoS One 14:e0216211. doi: 10.1371/journal.pone.0216211

PubMed Abstract | CrossRef Full Text | Google Scholar

DeSantis, T. Z., Hugenholtz, P., Larsen, N., Rojas, M., and Andersen, G. L. (2006). Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl. Environ. Microbiol. 72, 5069–5072. doi: 10.1128/aem.03006-05

PubMed Abstract | CrossRef Full Text | Google Scholar

Edgar, R. C. (2013). UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10:996. doi: 10.1038/NMETH.2604

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisler, M. C., Lee, M. R. F., Tarlton, J. F., Martin, G. B., Beddington, J., Dungait, J. A. J., et al. (2014). Agriculture: steps to sustainable livestock. Nature 507, 32–34.

Google Scholar

Favier, C. F., Vaughan, E. E., De Vos, W. M., and Akkermans, A. D. (2002). Molecular monitoring of succession of bacterial communities in human neonates. Appl. Environ. Microbiol. 68, 219–226. doi: 10.1128/AEM.68.1.219-226.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, H., Chi, X., Li, G., Qin, W., Song, P., Jiang, F., et al. (2020). Gut microbial diversity and stabilizing functions enhance the plateau adaptability of Tibetan wild ass (Equus kiang). Microbiol. Open 9, 1150–1161. doi: 10.1002/mbo3.1025

PubMed Abstract | CrossRef Full Text | Google Scholar

Georgios, O., Vieira, T. A. G., Carla, F., Lucas, B. M., Silva, M. V., Carvalho, B. R., et al. (2013). Fecal microbial diversity in pre-weaned dairy calves as described by pyrosequencing of metagenomic 16S rDNA. Associations of Faecalibacterium species with health and growth. PLoS One 8:e63157. doi: 10.1371/journal.pone.0063157

PubMed Abstract | CrossRef Full Text | Google Scholar

Gregory, C. J., Kyle, B., Bushman, F. D., DeSantis, Z., and Andersen, G. L. (2009). PyNAST: a flexible tool for aligning sequences to a template alignment. Bioinformatics 2, 266–267. doi: 10.1093/bioinformatics/btp636

PubMed Abstract | CrossRef Full Text | Google Scholar

Haas, B. J., Gevers, D., Earl, A. M., Feldgarden, M., Ward, D. V., Giannoukos, G., et al. (2011). Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR amplicons. Genome Res. 21, 494–504. doi: 10.1101/gr.112730.110

PubMed Abstract | CrossRef Full Text | Google Scholar

Harlow, B. E., Lawrence, L. M., and Flythe, M. D. (2013). Diarrhea-associated pathogens, Lactobacilli and cellulolytic bacteria in equine feces: responses to antibiotic challenge. Vet. Microbiol. 166, 225–232. doi: 10.1016/j.vetmic.2013.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Harlow, B. E., Lawrence, L. M., Hayes, S. H., Crum, A., and Flythe, M. D. (2016). Effect of dietary starch source and concentration on equine fecal microbiota. PLoS One 11:e0154037. doi: 10.1371/journal.pone.0154037

PubMed Abstract | CrossRef Full Text | Google Scholar

Heberle, H., Meirelles, G. V., Silva, F. R., Telles, G. P., and Minghim, R. (2015). InteractiVenn: a web-based tool for the analysis of sets through Venn diagrams. BMC Bioinform. 16:169. doi: 10.1186/s12859-015-0611-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Huber, T., Faulkner, G., and Hugenholtz, P. (2004). Bellerophon: a program to detect chimeric sequences in multiple sequence alignments. Bioinformatics 20, 2317–2319. doi: 10.1093/bioinformatics/bth226

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiao, J., Jian, W., Chuanshe, Z., Shaoxun, T., Min, W., and Zhiliang, T. (2016). Composition of Ileal bacterial community in grazing goats varies across non-rumination, transition and rumination stages of life. Front. Microbiol. 7:1364. doi: 10.3389/fmicb.2016.01364

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhl, J., Winterhoff, N., Wulf, M., Schweigert, F. J., Schwendenwein, I., Bruckmaier, R. M., et al. (2011). Changes in faecal bacteria and metabolic parameters in foals during the first six weeks of life. Vet. Microbiol. 151, 321–328. doi: 10.1016/j.vetmic.2011.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Langille, M. G. I., Zaneveld, J., Caporaso, J. G., McDonald, D., Knights, D., Reyes, J. A., et al. (2013). Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat. Biotechnol. 31, 814–821. doi: 10.1038/nbt.2676

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Xiaoxu, W., Ting, Z., Huazhe, S., Weixiao, N., Chao, X., et al. (2018). The development of microbiota and metabolome in small intestine of sika deer (Cervusnippon) from birth to weaning. Front. Microbiol. 9:4. doi: 10.3389/fmicb.2018.00004

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, G., Bou, G., Su, S., Xing, J., Qu, H., Zhang, X., et al. (2019). Microbial diversity within the digestive tract contents of Dezhou donkeys. PLoS One 14:e0226186. doi: 10.1371/journal.pone.0226186

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Fan, H., Ding, X., Hong, Z., Nei, Y., Liu, Z., et al. (2014). Analysis of the gut microbiota by high-throughput sequencing of the V5-V6 regions of the 16S rRNA gene in donkey. Curr. Microbiol. 68, 657–662. doi: 10.1007/s00284-014-0528-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Moreau, M. M., Eades, S. C., Reinemeyer, C. R., Fugaro, M. N., and Onishi, J. C. (2014). Illumina sequencing of the V4 hypervariable region 16S rRNA gene reveals extensive changes in bacterial communities in the cecum following carbohydrate oral infusion and development of early-stage acute laminitis in the horse. Vet. Microbiol. 168, 436–441. doi: 10.1016/j.vetmic.2013.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Mshelia, E. S., Adamu, L., Wakil, Y., Turaki, U. A., Gulani, I. A., and Musa, J. (2018). The association between gut microbiome, sex, age and body condition scores of horses in Maiduguri and its environs. Microb. Pathog. 118, 81–86. doi: 10.1016/j.micpath.2018.03.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Mueller, N. T., Bakacs, E., Combellick, J., Grigoryan, Z., and Dominguez-Bello, M. G. (2015). The infant microbiome development: mom matters. Trends Mol. Med. 21, 109–117. doi: 10.1016/j.molmed.2014.12.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Newell, P. D., and Douglas, A. E. (2014). Interspecies interactions determine the impact of the gut microbiota on nutrient allocation in Drosophila melanogaster. Appl. Environ. Microbiol. 80, 788–796. doi: 10.1128/AEM.02742-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Quercia, S., Freccero, F., Castagnetti, C., Soverini, M., Turroni, S., Biagi, E., et al. (2019). Early colonisation and temporal dynamics of the gut microbial ecosystem in Standardbred foals. Equine Vet. J. 51, 231–237. doi: 10.1111/evj.12983

PubMed Abstract | CrossRef Full Text | Google Scholar

Satokari, R., Grönroos, T., Laitinen, K., Salminen, S., and Isolauri, E. (2008). Bifidobacterium and Lactobacillus DNA in the human placenta. Lett. Appl. Microbiol. 48, 8–12. doi: 10.1111/j.1472-765X.2008.02475.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Schoster, A., Staempfli, H. R., Guardabassi, L. G., Jalali, M., and Weese, J. S. (2017). Comparison of the fecal bacterial microbiota of healthy and diarrheic foals at two and four weeks of life. BMC Vet. Res. 13:144. doi: 10.1186/s12917-017-1064-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Segata, N., Izard, J., Waldron, L., and Gevers, D. (2011). Metagenomic biomarker discovery and explanation. Genome Biol. 12:R60.

Google Scholar

Sonnenburg, E. D., Zheng, H., Joglekar, P., Higginbottom, S. K., Firbank, S. J., Bolam, D. N., et al. (2010). Specificity of polysaccharide use in intestinal Bacteroides species determines diet-induced microbiota alterations. Cell 141, 1241–1252. doi: 10.1016/j.cell.2010.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Stephens, W. Z., Burns, A. R., Stagaman, K., Wong, S., Rawls, J. F., Guillemin, K., et al. (2016). The composition of the zebrafish intestinal microbial community varies across development. ISME J. 10, 644–654. doi: 10.1038/ismej.2015.140

PubMed Abstract | CrossRef Full Text | Google Scholar

Torok, V. A., Ophelkeller, K., Loo, M., and Hughes, R. J. (2008). Application of methods for identifying broiler chicken gut bacterial species linked with increased energy metabolism. Appl. Environ. Microbiol. 74, 783–791. doi: 10.1128/aem.01384-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Uyeno, Y., Sekiguchi, Y., and Kamagata, Y. (2010). rRNA-based analysis to monitor succession of fecal bacterial communities in Holstein calves. Lett. Appl. Microbiol. 51, 570–577. doi: 10.1111/j.1472-765X.2010.02937.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Garrity, G. M., Tiedje, J. M., and Cole, J. R. (2007). Naive bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl. Environ. Microbiol. 73, 5261–5267. doi: 10.1128/aem.00062-07

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams, J. E., Carrothers, J. M., Lackey, K. A., Beatty, N. F., Brooker, S. L., Peterson, H. K., et al. (2019). Strong multivariate relations exist among milk, oral, and fecal microbiomes in mother-infant dyads during the first six months postpartum. J. Nutr. 149, 902–914. doi: 10.1093/jn/nxy299

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, G. D., Compher, C., Chen, E. Z., Smith, S. A., Shah, R., Bittinger, K., et al. (2016). Comparative metabolomics in vegans and omnivores reveal constraints on diet-dependent gut microbiota metabolite production. Gut 65, 63–72. doi: 10.1136/gutjnl-2014-308209

PubMed Abstract | CrossRef Full Text | Google Scholar

Yatsunenko, T., Rey, F. E., Manary, M. J., Trehan, I., Dominguez-Bello, M. G., Contreras, M., et al. (2012). Human gut microbiome viewed across age and geography. Nature 486, 222–227. doi: 10.1038/nature11053

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, W., Chow, K.-H., Fleming, E., and Oh, J. (2019). Selective colonization ability of human fecal microbes in different mouse gut environments. ISME J. 13:805. doi: 10.1038/s41396-018-0312-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Zoetendal, E. G., Akkermans, A. D. L., Vliet, W. M. A., De Visser, J. A. G. M., and DeVos, W. M. (2001). The host genotype affects the bacterial community in the human gastrointestinal tract. Microb. Ecol. Health Dis. 13, 129–134. doi: 10.1080/089106001750462669

CrossRef Full Text | Google Scholar