Comparative Analysis of the Gut Microbiota of Mongolian Gazelle (Procapra gutturosa) Under Fragmented Habitats

The existence of man-made facilities such as pasture fences makes the grassland ecosystem fragmented and endangers the survival of local wild animals. The Mongolian gazelle is highly sensitive to hunting and habitat destruction, and is one of the most threatened artiodactyls in Eurasia. It provides a critical model to studying gut microbiota under fragmented habitats. Therefore, we applied metagenomics sequencing to analyze the gut microbiota communities and functions of Mongolian gazelle under fragmented habitats. The results demonstrated that there were no significant differences in gut microbial communities between the different groups at both the phylum and genus level. The functional analyses showed that the Mongolian gazelle in fragmented habitat had a stronger ability to degrade naphthalene, but their ability to absorb carbohydrates was weaker. This study provided fundamental information about the gut microbiota of Mongolian gazelle, and we recommend reducing habitat fragmentation to better protect the Mongolian gazelle.


INTRODUCTION
The gut microbiota of mammals is frequently affected by many factors, including dietary choices, phylogeny, and environmental changes (Turnbaugh et al., 2009;Yatsunenko et al., 2012;Schnorr et al., 2014;Clemente et al., 2015). The diverse and extremely complex gut microbiota benefits host in many ways, such as synthesizing vitamins, stimulating immune responses, and performing metabolic functions that the host cannot perform (Backhed et al., 2005;Cerf-Bensussan and Gaboriau-Routhiau, 2010;Leblanc et al., 2013;Crespo-Piazuelo et al., 2018;Nagarajan et al., 2018). Ruminants rely on the gut microbiota to harness energy via the fermentation of dietary material (Thoetkiattikul et al., 2013). Thus, the gut microbiota plays an important role in the nutritional ecology of ruminants . Recent studies indicate that anthropogenic disturbance can cause major changes in the composition of the gut microbiota (Barelli et al., 2020). For instance, populations of both the black howler monkey and Udzungwa red colobus have been shown to have lower gut microbiota diversity in fragmented habitats compared with intact habitats (Amato et al., 2013;Barelli et al., 2015). Given that the gut microbiota plays a crucial role in host health (Fung et al., 2017), understanding the effects of anthropogenic disturbance on host's gut microbiota is critical for developing effective conservation strategies for endangered species (Stumpf et al., 2016).
The Mongolian gazelle (Procapra gutturosa) is a unique herbivorous animal on the steppe habitats in Eurasia and was listed as a Category I species in the China's Red List of Biodiversity: Vertebrates (Zhang et al., 2014;Jiang et al., 2021). The species was once widespread in Mongolia, northern China, and southeastern Siberia (Wang et al., 1997). A previous study of Mongolian gazelle has focused on behavioral characteristics, feeding habits, and migration (Li et al., 1999;Leimgruber et al., 2001;Odonkhuu et al., 2009). Mongolian gazelle population in Hulun Lake National Nature Reserve is affected by anthropogenic disturbance due to the existence of human facilities such as grassland fence, and the habitat is fragmented. Affected by the fragmented environment, the Mongolian gazelle formed a locally isolated population in this area, which greatly increased the risk of its local extinction. In contrast, the Mongolian gazelle population in the China-Mongolia border area has been virtually unaffected by anthropogenic disturbance, and the habitat is complete. Therefore, we speculate that the composition and function of gut microbiota would be altered in populations restricted to fragmented habitats.
The present study aimed to explore gut microbial diversity and function of Mongolian gazelle populations under fragmented habitats by using the metagenomic sequencing. This research will improve the understanding of the differences in gut microbiota of the Mongolian gazelle in different conditions and provide a scientific reference for the protection and management of this species.

Sample Collection
Samples were collected from China-Mongolia border area and Hulun Lake National Nature Reserve. Detailed information for all samples is shown in Supplementary Table 1. Eight  microsatellite loci (OArFCB304, SPS115, TGLA68, IOBT395,  PZE114, MNS72, BM1341, and MB066) were used to identify individuals (Buchanan and Crawford, 1993;Reed and Beattie, 2001;Chu et al., 2002;Motavalian et al., 2002;Lv et al., 2008;Fend et al., 2010;Nijman et al., 2010). This identification allows all alleles identical or only one mismatch (Ning et al., 2018). It was identified that the fecal samples of H group (n = 4) and B group (n = 5) were all from different individuals. Detailed information for microsatellite loci of Mongolian gazelle is shown in Table 1. During the sampling period, the ambient temperature was approximately −30 • C to ensure the quality of DNA from gut microbiota. Fecal samples were placed in sterile containers and stored at −80 • C until DNA extraction.

DNA Extraction, Library Construction, and Metagenomics Sequencing
According to the manufacturer's recommendations, metagenomic DNA was extracted from the internal part of feces using QIAamp Fast DNA Stool Mini Kit (Qiagen, Germany).
Extracted DNA was monitored using 1% agarose gels to determine degradation degree and potential contamination. We measured the concentration of DNA with Qubit dsDNA Assay Kit in Qubit 2.0 Flurometer (Life Technologies, CA, United States). Only samples that meet the following criteria were used to construct the library: (1) OD value is between 1.8 and 2.0; (2) DNA contents above 1 µg.
The qualified DNA was fragmented to a size of 350 bp using sonication. DNA fragments were end-polished, A-tailed, and ligated with the full-length adaptor for sequencing and further PCR amplification. The PCR products were purified with AMPure XP system. Next, we analyzed the size distribution of libraries by Agilent 2100 Bioanalyzer and quantified libraries by real-time PCR. According to the manufacturer's instructions, we performed the clustering of the index-coded samples on a cBot Cluster Generation System. After cluster generation, Illumina HiSeq platform was used to sequence the library preparations and paired-end reads were generated.
We blasted the unigenes to the sequences of bacteria, fungi, archaea, and viruses from the NR database of NCBI

Summary of the Sequencing Data
A total of 114 Gb raw data were obtained from the samples. To ensure the accuracy and reliability of downstream analysis, we carried out quality control and filtering. The total amount of clean data also remained at about 114 Gb, and the average efficiency was more than 99.8% ( Table 2).

Gene Prediction and Abundance Analysis
To evaluate whether the collected samples could meet the requirements of subsequent bioinformatics analysis, we conducted a rarefaction curve analysis. The results of rarefaction curve analysis based on the core genes ( Figure 1A) and pan genes ( Figure 1B) showed that our gene catalog captured all the available gene information in our sample. The heatmap of correlation coefficients showed that the correlation between groups is smaller than the correlation within groups, indicating reliable experimentation and reasonable sample selection ( Figure 1C). From the box-plot diagram, we found that the number of genes was higher in B group compared with that in H group ( Figure 1D).
Human exploitation and destruction of resources are currently threatening innumerable wild animal species, altering natural ecosystems and, thus, food resources, with profound effects on gut microbiota (Barelli et al., 2020). To assess the differences in gut microbiota as affected by the anthropogenic disturbance, we applied the Metastats analysis. H group and B group showed significant differences (q < 0.05) in the following species: Bacteroides sp. HPS0048, Chryseobacterium sp. CBo1, Afipia broomeae, Pseudoscardovia radai, Microcystis phage MaMV-DC, Chloroflexi bacterium CG_4_9_14_3_um_filter_45_9, Anaerobacillus macyae, and Euryarchaeota archaeon SM23-78 (Figure 3). PCA at the species level showed that H group and B group have distinct cluster regions (ANOSIM: p = 0.014, R = 0.444) (Supplementary Figure 1). It is remarkable that the gut microbiota did not differ significantly between H group and B group at the both phylum and genus level. Anthropogenic disturbance only caused some changes in the gut microbiota of Mongolian gazelle at the species level.

Common Functional Database Annotations
To explore the activity of genes in gut microbiota, we annotated the functional and metabolic pathways based on the KEGG and CAZy databases (Figure 4).
Based on the KEGG database, we found that the gut microbiota of Mongolian gazelle has enriched activity for metabolism of carbohydrates (gene number: 92205), amino acid (68527), nucleotide (48992), and cofactors and vitamins (42745) (Figure 4A). The third level of KEGG classification indicated that the relative abundance of naphthalene degradation (ko00626) was significantly higher in H group than in B group (q < 0.05) (Figure 5A). Through Metastats analysis, we found that the alcohol dehydrogenase (EC:1.1.1.1) was significantly higher in H group than in B group (q < 0.05) ( Figure 5B). The enzyme is mainly involved in biological processes, such as fatty acid degradation; glycine, serine, and threonine metabolism; microbial metabolism in diverse environments; and biosynthesis of secondary metabolites. For the H group, alcohol dehydrogenase was related to naphthalene degradation (Supplementary Figure 2).
From the first level of CAZy classification, we found that most of the genes were annotated to the three functional configurations: glycoside hydrolases (GHs), glycosyltransferases (GTs), and carbohydrate-binding modules (CBMs) (Figure 4B). At the second classification level, GH2 (B group: 0.28%; H group: 0.24%), GT2 (B group: 0.23%; H group: 0.24%), and GH43  (Supplementary Figure 3). Among families, 6 were significantly more abundant in H group than in B group (q < 0.05), including GH65/121/24/64/113 and CBM56; 7 were significantly more abundant in B group than in H group (q < 0.05), including GH138/142/106/2, GT5, polysaccharide lyase family 1 (PL1), and CBM62 (Figure 6). At the third classification level, a total of 35 enzymes showed significant differences (q < 0.05) (Supplementary Figure 4). Twenty-eight enzymes were significantly more abundant in B group than in H group, including 22 enzymes of GH, 5 enzymes of GT, and 1 enzyme of PL; 7 enzymes of GH were significantly more abundant in H group than in B group.

DISCUSSION
Complete habitats offer wild animal inhabitants a more diverse diet than fragmented habitats. Furthermore, it has been suggested that dietary variation promotes the changes in gut microbiota (Clarke et al., 2014;Corlett, 2016). In this study, we compared the differences in gut microbiota of the Mongolian gazelle in the different habitat types. Our research showed that Firmicutes and Bacteroidetes were the prominent phyla in the Mongolian gazelle, which is consistent with previous studies of other ruminants (Sundset et al., 2007;Guan et al., 2017;Li et al., 2017). Firmicutes can degrade fibers into short-chain fatty acids, while the main functions of Bacteroidetes are to degrade fats, proteins, and carbohydrates (Jami et al., 2013). Ruminants rely on gut microbiota to harness energy by fermenting dietary material (Thoetkiattikul et al., 2013). Hence, the higher abundance of Firmicutes and Bacteroidetes can be correlated with the dietary choices of Mongolian gazelle. At the species level, we found that Bacteroides sp. HPS0048 was significantly more abundant in B group than in H group. The higher abundance of Bacteroides might indicate a healthier gut microecosystem because it can promote the improvement of immune system and maintain the balance of gut microbiota (Hooper et al., 2001;Hooper, 2004;Sears, 2005). We observe that the habitat fragmentation had no significant effect on the gut microbiota of Mongolian gazelle at both the phylum and genus level. The core gut microbiota of Mongolian gazelle at the different habitat types was relatively stable, suggesting that the effects of habitat fragmentation on gut microbiota were limited. Previous studies on ruminants indicated that 32 species had relatively stable core gut microbiota regardless of variations in dietary choices (Henderson et al., 2016). These results all indicated that ruminants had relatively stable core gut microbiota, and environment might affect the abundance of gut microbiota (Weimer, 2015).
To further understand the environmental adaptations of Mongolian gazelle, we carried out the functional analyses of gut microbiota. From the KEGG analyses, it was obvious that the metabolic pathway of naphthalene degradation has higher abundance in H group. Naphthalene belongs to polycyclic aromatic hydrocarbons and is one of the most prevalent contaminants (Mihelcic and Luthy, 1991;Ghoshal et al., 1996;Marx and Aitken, 1999). Previous studies have shown that the gut microbiota of snails can degrade naphthalene (Hu Z. et al., 2018). Given that H group is    more affected by anthropogenic disturbance, we speculate that the higher abundance of naphthalene degradation is essential to maintain the stability of gut microbiota in the Mongolian gazelle.
Based on the analysis of CAZy database, we found that the GH families and CBM families have larger proportion differences between the two populations of Mongolian gazelle. The relative abundance of GH65/121/24/64/113 and CBM56 were higher in H group, whereas GH138/142/106/2 and CBM62 were higher in B group. GH families can facilitate the hydrolysis of cellulose so that Mongolian gazelle had stronger capacity for fiber digestion (Hess et al., 2011). CBM families do not show enzymatic activity, but can help GH families bind to polysaccharides and strengthen their activity (Du et al., 2018;Bernardes et al., 2019). GH families and CBM families of Mongolian gazelle can help them maximize energy from the cellulose, ensuring them to adapt to the wild environments. Our data also showed that PL1 and GT5 were significantly more abundant in B group than in H group. PL1 can cleave glycosidic bonds of homogalacturonan (HG), which is a multifunctional pectic polysaccharide of plant cell walls, and GT families can use sugar donors containing nucleoside phosphate or lipid phosphate leaving groups to catalyze the formation of glycosidic bonds (Willats et al., 2001;Lairson et al., 2008;Abbott et al., 2013). At the third classification level, we also observed that five enzymes belonging to GT5 (UDP-Glc: alpha-1,4-glucan synthase, ADP-Glc: starch glucosyltransferase, UDP-Glc: alpha-1,3-glucan synthase, UDP-Glc: glycogen glucosyltransferase, and NDP-Glc: starch glucosyltransferase) were significantly more abundant in B group than in H group. Therefore, the high abundance of PL1 and GT5 in B group showed their stronger ability to absorb carbohydrates.

CONCLUSION
This study characterized the gut microbiota of Mongolian gazelle under fragmented habitats using metagenomics sequencing. Compositions of gut microbiota were similar between H group and B group, but the functions of gut microbiota differed. In the complete habitat, the high relative abundance of PL1 and GT5 promoted the absorption of carbohydrates. The higher abundance of naphthalene degradation in H group could enable Mongolian gazelle to maintain the stability of gut microbiota in the habitat that was more affected by anthropogenic disturbance. After comparison, we speculated that the variation in functions of gut microbiota in fragmented habitat might be from anthropogenic disturbance. Therefore, we recommended that minimal human disturbance and complete habitats were crucial for the survival of Mongolian gazelle.

ETHICS STATEMENT
The animal study was reviewed and approved by the Qufu Normal University Institutional Animal Care and Use Committee.

AUTHOR CONTRIBUTIONS
HZ, XY, and LS conceived and designed the study. LS, XY, HD, TL, LW, SZ, YS, and YD performed the research. LS and XY analyzed the data and prepared the manuscript. All authors read and approved the final article.

FUNDING
This work was supported by the National Natural Science Foundation of China (31872242, 32070405, 31900311, 32000291, 32001228, and 32170530).

ACKNOWLEDGMENTS
We thank staff of Hulun Lake National Nature Reserve and World Wide Fund for Nature for helping with sample collection.