Gut Microbiome Changes in Captive Plateau Zokors (Eospalax baileyi).

Wild-caught animals must cope with drastic lifestyle and dietary changes after being induced to captivity. How the gut microbiome structure of these animals will change in response receives increasing attention. The plateau zokor (Eospalax baileyi), a typic subterranean rodent endemic to the Qinghai-Tibet plateau, spends almost the whole life underground and is well adapted to the environmental pressures of both plateau and underground. However, how the gut microbiome of the plateau zokor will change in response to captivity has not been reported to date. This study compared the microbial community structure and functions of 22 plateau zokors before (the WS group) and after being kept in captivity for 15 days (the LS group, fed on carrots) using the 16S rRNA gene via high-throughput sequencing technology. The results showed that the LS group retained 973 of the 977 operational taxonomic units (OTUs) in the WS group, and no new OTUs were found in the LS group. The dominant bacterial phyla were Bacteroides and Firmicutes in both groups. In alpha diversity analysis, the Shannon, Sobs, and ACE indexes of the LS group were significantly lower than those of the WS group. A remarkable difference (P < 0.01) between groups was also detected in beta diversity analysis. The UPGMA clustering, NMDS, PCoA, and Anosim results all showed that the intergroup difference was significantly greater than the intragroup difference. And compared with the WS group, the intragroup difference of the gut microbiota in the LS group was much larger, which failed to support the assumption that similar diets should drive convergence of gut microbial communities. PICRUSt revealed that although some functional categories displayed significant differences between groups, the relative abundances of these categories were very close in both groups. Based on all the results, we conclude that as plateau zokors enter captivity for a short time, although the relative abundances of different gut microbiota categories shifted significantly, they can maintain almost all the OTUs and the functions of the gut microbiota in the wild. So, the use of wild-caught plateau zokors in gut microbial studies is acceptable if the time in captivity is short.

Introduction

Plateau zokors (Eospalax baileyi) belong to Spalacidae, Rodentia,[1] and are a typical species of subterranean herbivores. They are endemic to the Qinghai-Tibet plateau (QTP) and the surrounding high-altitude areas, and mainly dwell in alpine meadows, alpine grasslands, shrubland, and farmland in the eastern QTP.[2] As food competitors with livestock and agents of soil erosion in QTP, plateau zokors have traditionally been viewed as pests.[3] After a long time survival in such a habitat, the plateau zokor is well adapted to the environmental pressures induced by high altitude, such as solar radiation,[4] hypoxia,[5,6] and low temperatures.[7-10] As a typical subterranean rodent, the plateau zokor lives underground for almost its entire life,[11] where most activities (including foraging and reproduction) are carried out in the subterranean burrow. Because of the limitations imposed by underground burrows and the influence of seasonal plant withering, in winter,[12-15] the plateau zokor is faced with a shortage of food resources and fewer chances to choose what to eat,[16] which should be responsible for its habit as a dietary generalist[17] and preference even for some common poisonous weeds, such as Oxytropis kansuensis[18] and Stellera chamaejasme.[17] Benefiting from all the characteristics described above, the plateau zokor is an ideal model to study hypoxia adaptation, low temperature adaptation, and detoxification mechanisms.[19,20] As a hindgut fermenter, the plateau zokor is supposed to be closely associated with bacteria.[21] However, as far as we know, there is little research on its symbiotic microbiota.

Numerous studies have shown that many animals are associated with symbiotic communities of microorganisms.[22] These microbiotas are closely associated with nutrition, metabolism, immunity, and health of the host.[23-25] In fact, the microbiome has been referred to as the “second genome” of animals.[26-30] With regard to the great effect the gut microbiota exerts on the host, studying the gut microbiota of plateau zokor provides new insight into its special adaptive mechanism. Additionally, the development of PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) greatly facilitates the prediction of bacterial functions based on 16S rRNA gene sequencing data.[31] However, the diversity of the host gut microbiota can be affected by a number of different factors, such as species, diet, genotype, geography, and age.[32-36] Recent studies increasingly reported the effects of captivity on the microbiome of wild animals.[33,37-39] Since these animals are exposed to a series of changes in lifestyle, after being transferred into an artificial environment from the wild, the gut microbiota of wild animals will change greatly.[40,41] However, the trends of these changes are not always identical among species, and little is known about the critical factors that are responsible for these changes.[40] Thus, far more comparative studies of gut microbiota between captive and wild members of a species are needed to fill some important knowledge gaps.

Therefore, it is important to gain a perspective on the differences of the gut microbiota of the plateau zokor in captivity versus their wild state and incorporate this knowledge into future animal microbiome study design. In the case of the plateau zokors, after their capture, many researchers will feed them with carrots in the laboratory for about 2 weeks,[42] so that the plateau zokor can adapt to the new environment. However, after being captured, plateau zokors were transferred from their underground burrows to the ground, and thus both their living environment and diet undergo a drastic change. The changes that happen to the gut microbiota of plateau zokor in response remain unclear. However, until now, the gut microbiota of the plateau zokor has not been reported extensively, with the exception being the research on the isolation of tannin-degrading bacteria from the cecum of the plateau zokor.[43] This study compared the changes of the gut microbiota and bacterial functions of 22 male plateau zokors in response to feeding in the laboratory with captive diet with their wild state, utilizing the 16S rRNA gene and high-throughput sequencing technology.

Furthermore, food is a key factor to determine the gut microbiota composition.[44-46] Similar diets appear to drive convergence of gut microbial communities.[37,47] Therefore, while designing control experiments on the gut microbiota, according to the requirement of a single variable principle, researchers commonly choose to provide the experimental subjects with the same diet for a specific time after bringing them into the laboratory; consequently, their gut microbiota composition will become stable and consistent before the beginning of control experiments.[48,49] In view of this, the present study also evaluated whether the microbiota of captive plateau zokors will converge after individuals are fed a diet based on roots (carrots).

Materials and Methods

Animals and diet

Twenty-two male plateau zokors were captured with live traps in May 2019 in Datong county (37°8′20″N, 101°15′1″E, elevation 3111 m), Qinghai, China. The animals were maintained individually in stainless steel cages (40 cm × 30 cm × 25 cm, labeled 1-22) in a constant environment (at 25°C ± 1°C, wood shavings were used as bedding material, and the windows were shaded with shade curtains). To ensure the same diet, all the plateau zokors were fed with carrots only as the captivity diet. According to our previous experiment, plateau zokors raised with carrots only can live well, and each zokor consumes about 150 g carrot a day. Here, to ensure enough food, each plateau zokor was provided with 250 g fresh carrot daily and the remains of the day before were removed.

Sample collection

As soon as the plateau zokors were brought to the laboratory, fresh fecal samples were initially collected from each plateau zokor, labeled W01 to W22 (ie, the WS group). Then, the plateau zokors were fed with carrots for 15 days, and fresh fecal samples were collected from each plateau zokor on day 15, labeled L01 to L22 (ie, the LS group). In addition, all samples were frozen in liquid nitrogen as soon as possible after collection. All procedures involving the handling and care of animals were in accordance with the China practice for the care and use of laboratory animals. The experimental design and procedures were approved by the Animal Care and Use Committee of Northwest Institute of Plateau Biology, Chinese Academy of Sciences (IACUC Issue No. NWIPB2019009).

DNA extraction and sequencing

DNA from all samples was extracted using the cetyl trimethyl ammonium bromide (CTAB) method.[50] To obtain 16S V3 and V4 regions, target genes were amplified, and using a specific barcode primer set (341F, 806R).[51] All PCR reactions were performed in a 30 μL reaction mixture (15 μL Phusion Master Mix 2X, New England Biolabs; 3 μL [6 μM] primer [2 μM]; 10 μL [5-10 ng] DNA [1 ng/μL]; 2 μL dd H2O). Thermal cycling consisted of an initial denaturation at 98°C for 1 minute, followed by 30 cycles of denaturation at 98°C for 10 seconds, annealing at 50°C for 30 seconds, and extension at 72°C for 30 seconds; a final extension was performed at the end of the run at 72°C for 5 minutes. An identical volume of 1X loading buffer was mixed with the PCR products; then, electrophoresis was performed on a 2% agarose gel for PCR product detection. Bright bands (400-450 bp) were chosen for further experiments. PCR products were mixed at equidensity ratios, and were then extracted with the GeneJET Gel Extraction Kit (Thermo Scientific).

Library preparation and sequencing

The Illumina TruSeq DNA PCR-Free Library Preparation Kit (Illumina, USA) was used to generate sequence libraries according to the manufacturer’s recommendations and index codes were added. A Qubit@ 2.0 Fluorometer (Thermo Scientific) and an Agilent 2100 Bioanalyzer system were used to assess the quality of the generated libraries. Once they were qualified, the library was sequenced on an Illumina NovaSeq 6000 platform.

Data analysis

Single-end reads were assigned to samples based on their unique barcodes and were truncated by cutting off both barcodes and primer sequences. Quality filtering on the raw reads was performed under specific filtering conditions: (1) The 300 bp reads were truncated at any site that received an average quality score of <20 over a 50 bp sliding window. Truncated reads shorter than 50 bp or containing ambiguous characters were discarded. (2) Only overlapping sequences longer than 10 bp were assembled according to their overlapped sequence. The maximum mismatch ratio of the overlap region was 0.2. Reads that could not be assembled were discarded. (3) Samples were distinguished according to their barcode and primers, and the sequence direction was adjusted. Then, high-quality clean reads were obtained. These reads were compared with the reference database using the UCHIME algorithm to detect and then remove chimera sequences.[52,53] Sequence analyses were conducted by Uparse software.[54] Sequences with a similarity ⩾97% were assigned to a same operational taxonomic unit (OTU). OTUs that were identified in the blank controls were removed. OTUs that only appeared in 1 sample and OTUs with a total sequence count of <5 were also removed. The representative sequence for each OTU was screened for further annotation. The taxonomy of each OTU representative sequence was analyzed by Mothur[55] against the SILVA 132/16s_bacteria database[56] with the confidence threshold of 0.8.

The OTU table was rarified and 4 alpha diversity metrics were calculated: Sobs (observed richness) index, Shannon index, ACE, and Goods coverage. Rarefaction curves were generated based on these 4 metrics. QIIME[57] was used to calculate both weighted and unweighted Unifrac,[58] which are phylogenetic measures of beta diversity. Unweighted unifrac was used for Unweighted Pair Group Method with Arithmetic (UPGMA) mean clustering. LEfSe[59] analysis was conducted by LEfSe 1.0 with a defaulted filter value of 4 for the LDA score. Principal coordinates analysis (PCoA) was displayed using the vegan package, stat package, and ggplot2 package.[60] For non-metric multi-dimensional scaling (NMDS) analysis, the vegan package for R[61] was used. STAMP[62] was used for comparing the differences between groups at different categories. Bacterial functions were predicted based on 16S rRNA gene by the PICRUSt algorithm.[31] Considering the limitations of PICRUSt[63] and its use on the non-model organism, we suggest some cautions on the interpretation of the results.

Novogene Co., Ltd., was commissioned to complete all experiments described in the sections on DNA extraction and PCR amplification, library preparation and sequencing, and previous data analysis.

Results

Body weight

After feeding the Plateau zokors with carrots for 15 days, the average body weight changed from 325.45 ± 65.21 g (WS) to 335.18 ± 57.83 g (LS). Paired t-tests indicated a significant increase in body weight (t = −2.834, P = 0.010).

Data quality

A total of 3 921 111 reads of raw data were obtained, and the average number of reads per sample was 89 116 ± 6389. The combined reads were 3 643 113, accounting for 92.91%. A total of 3 240 868 qualified reads with an average length of 418.53 bp were included in the clean data, with an average effective ratio of 88.96%. The Q20 percentage was 97.96% and the Q30 percentage was 93.78%. Based on the qualified reads, 977 OTUs were harvested, 973 of which were detected in both the wild and captive groups while the remaining 4 and 0 OTUs were detected only in the WS or LS groups, respectively. The rarefaction curves of sobs (the quantity of observed OTUs) and Shannon index values on OTUs level became gradually placid as the sequencing depth increased. Both indexes demonstrated that each fecal sample of the plateau zokor had sufficient OTUs to reflect the maximum level of bacterial diversity, and almost all the bacterial species were detected at present sequencing depth. Judging from the gradual horizontalization of the rarefaction curve and rank abundance (see Supplemental Figure 1), the sample size and sequencing depth were adequate.

Gut microbiota profile

Other than unclassified, in total, all identified OTUs were sorted into 11 phyla, 23 orders, 41 families, and 97 genera. And comparison of the numbers of different microbiological taxonomic units of WS and LS groups is shown in Table 1.

Table 1.

Numbers of different microbiological taxonomic units in this study.

Group	Phylum	Class	Order	Family	Genus	OTU
LS	11	16	23	41	97	973
WS	11	16	23	41	97	977
Total	11	16	23	41	97	977

The relative abundances at the phylum level are shown in Figure 1. Bacteroidetes and Firmicutes were predominant phyla in both WS and LS groups, with Bacteroidetes counting for 58.49% in the WS group and 75.82% in the LS group, which were far more abundant than Firmicutes (38.36% in the WS group and 22.62% in the LS group). According to the Wilcoxon rank-sum test, the relative abundances of 10 phyla differed significantly between groups as shown in Supplemental Figure 2.

Figure. 1.

Relative abundance histograms: (a) The relative abundances for WS (wild group) and LS (captive group), and (b) The relative abundances for all samples. Others represents the phyla with relative abundance less than 1%.

At the family level, Wilcoxon rank-sum test indicated that 33 families were significantly different between both groups (see Supplemental Figure 3). After being fed in the laboratory with carrots, the relative abundance of 28 families in the LS group were significantly lower than those in the WS group, while 5 families were significantly higher. Furthermore, the relative abundances of the top 4 families in the LS group were Muribaculaceae (75.79% ± 12.73%), Lactobacillaceae (9.00% ± 13.13%), Ruminococcaceae (8.35% ± 5.21%), and Lachnospiraceae (3.57% ± 2.25%); while the abundances of these 4 families in the WS group were 58.39% ± 10.94% (Muribaculaceae), 0.49% ± 1.24% (Lacto-bacillaceae), 19.26% ± 5.76% (Ruminococcaceae), and 16.57% ± 6.02% (Lachnospiraceae). Compared with the WS group, in the LS group, Muribaculaceae increased by 29.80%, Lactobacillaceae increased by 1737%, Ruminococcaceae decreased by 56.65%, and Lachnospiraceae decreased by 78.46%.

At the genus level, the relative abundances of 81 genera were significantly different between groups (see Supplemental Figure 4).

Analysis of intragroup and intergroup differences

At the alpha diversity, the Goods coverages of the WS group and the LS group were both 99.8%, indicating a high level of diversity coverage. The Shannon index, ACE, Sobs, and Good coverage for both groups are shown in Table 2 and Figure 2. All these indexes supported that the diversity of the gut microbiota communities of the WS group was higher than that of the LS group.

Table 2.

Alpha diversity comparison between both groups based on multiple indexes.

Estimators	LS		WS		P-value	Q-value
	Mean	SD	Mean	SD
Shannon	4.152	0.462	5.076	0.168	0.000	0.000
ACE	749.709	74.388	889.990	25.775	0.000	0.000
Sobs	670.632	65.952	802.000	101.636	0.000	0.000
Goods coverage	0.998	0.000	0.998	0.000	0.000	0.000

Figure. 2.

Diversity boxplots with points for observed richness: (a) and Shannon diversity and (b) WS represents the wild group, LS represents that captive group.

As shown in the UPGMA clustering tree (Figure 3a), the NMDS (Figure 3c) and PCoA (Figure 3d) plots, based on the composition of the gut microbiota communities, samples from WS and LS groups were apparently separated, while the difference among individuals of the LS group were larger than those of the WS group. In the result of Anosim analysis (Figure 3b), the intergroup difference was significantly larger than the intragroup difference and the intragroup difference of the LS group was larger than that of the WS group (R = 0.778, P = 0.001). This is in complete agreement with the results of the NMDS and the UPGMA clustering analyses.

Figure. 3.

Beta diversity analysis results: (a) unweighted pair-group method with arithmetic mean (UPGMA) tree of unweighted unifrac distances. (b) Analysis of similarity (ANOSIM) between groups, (c) Bray_curtis distance-based non-metric multidimensional scaling (NMDS) analysis, and (d) Principal coordinate analysis of an unweighted unifrac distance matrix (PCoA). WS represents the wild group, LS represents that captive group.

In the LEFSe analysis (Figure 4), 16 microbial taxa were significantly different between groups, with 2 at the phylum level (Firmicutes and Bacteroidates), 3 at the class level (Clostridia, Bacilli, and Bacteroidia), 3 at the order level (Clostridiales, Lactobacillales, and Bacteroidales), 4 at the family level (Lachnospiraceae, Ruminococcaceae, Lactobacillaceae, and Muribaculaceae), and 4 at the genus level (unclassified_f_Lachnospiraceae, Ruminiclostridium_6, Lactobacillus, and norank_f_Muribaculaceae).

Figure 4.

Results of LDA effect size (LEfSe) analysis. Histogram of the LEfSe bar: (a) showed significant differences among the microbial taxa between WS and LS groups. The length of the column (ie, the LDA score) represents the influencing degree of each taxon. In the cladogram, and (b), the circle radiating inside-out represents the classification of the phylum to the species level. WS represents the wild group, LS represents that captive group.

Prediction of the functional composition

We predicted functions of gut bacterial communities using PICRUSt based on Kyoto Encyclopedia of Genes and Genomes (KEGG), http://www.genome.jp/kegg/. According to the level 1 KEGG modules (Table 3), the relative abundance of predicted functions in both groups varied in the same order: Metabolism > Genetic Information Processing > Environmental Information Processing > Cellular Processes > Organismal Systems > Human Diseases. And the relative abundance values of all the functional categories in the LS group were very close to the WS group, although the relative abundance of the other 5 functional categories, except for the Organismal Systems category, were detected to be significantly different between groups. According to the level 2 KEGG modules (see Supplemental Table 1), of all 39 detected functional categories, the relative abundance of 22 categories were higher than 1%. The relative abundances of these categories were also very close in both groups although some categories displayed significant differences between groups.

Table 3.

Relative abundance of predicted function for specific KEGG modules (level 1).

KEGG modules level 1	Relative abundance for specific KEGG modules (%)				P-value
	WS		LS
	Mean	SD	Mean	SD
Cellular processes	3.630	0.243	3.110	0.269	0.000
Environmental information processing	14.500	0.517	13.710	0.610	0.000
Genetic information processing	20.690	0.308	21.400	0.549	0.000
Human diseases	0.700	0.032	0.750	0.083	0.012
Metabolism	46.130	0.433	46.860	0.499	0.000
Organismal systems	0.780	0.030	0.780	0.040	0.881
Unclassified	13.570	0.140	13.390	0.301	0.015

Discussion

This study investigated the gut microbiota of the plateau zokor for the first time and compared the changes of gut microbiota before and after laboratory feeding. At the phylum level, the predominant components in the gut microbiota of both WS and LS groups were Bacteroidetes and Firmicutes. This is consistent with a large number of studies on herbivores,[64-66] indicating that the plateau zokor exhibits a high degree of adaptation to herbivores. In the LS group, totally 973 OTUs were detected, all of which were retained from the 977 OTUs in the WS group and accounted for 99.59% of the latter. It means that after being transferred to artificial environment, plateau zokors received very few microbes from food and environment, which is consistent with the study on woodrats.[67] And, it was reported that compared to the dietary specialists, food generalists are more likely to maintain the stability of gut microbiota,[33] just as the case of our study here on plateau zokors, a food generalist confirmed by previous studies.[17] This indicated that after plateau zokors enter captivity for a short time, they can retain almost all the OTUs of the gut microbiota in the wild. However, according to previous studies, compared with wild state, the gut microbiota of wild animals will change significantly in captivity because of the changes in lifestyle and diet.[39,40] Here, the plateau zokor is no exception. In the alpha diversity analysis, Shannon, Sobs, and ACE indexes of the WS group were significantly higher than those of the LS group. And a remarkable difference (P < 0.01) between groups was also detected in beta diversity analysis. Since the 2 groups shared almost all their OTUs, we think the detected difference between the 2 groups may have resulted from a shift in the relative abundances of different gut microbiota communities. Meanwhile, as shown in the function prediction by PICRUSt, according to both the KEGG modules level 1 and level 2, although some functional categories displayed significant differences between groups, the relative abundances of these categories, as well as all other categories detected, were very close in both groups. Based on the similarity in the OTU composition and the function of gut microbiota between the LS and WS groups, we conclude that the use of wild-caught plateau zokors for gut microbial studies is acceptable, at least for short periods of captivity.

Food has been frequently reported to be the main factor affecting the gut microbiota community of the host.[68-70] At the same time, many researchers also assume that the same diet can lead to the convergence of the gut microbiota of different individuals (both interspecific and intraspecific).[37,46,47] As the plateau zokors were subjected into captivity, they have to cope with drastic lifestyle and dietary changes. If diet plays a key role in the formation of microbial communities, then it would be expected that the microbial communities of plateau zokors would converge after exposure to a same diet. However, the UPGMA clustering, NMDS, PCoA, and Anosim results showed that, compared with the LS group, the gut microbial compositions of different plateau zokor individuals in the WS group were much more similar, which failed to support the assumption that similar diets should drive convergence of gut microbial communities. Consequently, we think that diet may not be the dominant factor affecting the gut microbial compositions of plateau zokors after being captured, but more further research is needed to address this question.

It has been reported that the plateau zokor prefers some common poisonous weeds, such as Oxytropis kansuensis[18] and Stellera chamaejasme,[17] which can cause intoxication in livestock and thus do great harm to pasture animal husbandry. Many researchers believe that the gut microbiota can contribute to the degradation of plant secondary metabolites for hosts.[27,30,71] So, further study on the gut microbiota of the plateau zokor may help us to understand the detoxification mechanism of the plateau zokor, which may provide new insight into the solution of poisonous weeds problem.

Click here for additional data file.

Supplemental material, sj-tif-1-evb-10.1177_1176934321996353 for Gut Microbiome Changes in Captive Plateau Zokors (Eospalax baileyi) by Daoxin Liu, Pengfei Song, Jingyan Yan, Haijing Wang, Zhenyuan Cai, Jiuxiang Xie and Tongzuo Zhang in Evolutionary Bioinformatics

Click here for additional data file.

Supplemental material, sj-tif-2-evb-10.1177_1176934321996353 for Gut Microbiome Changes in Captive Plateau Zokors (Eospalax baileyi) by Daoxin Liu, Pengfei Song, Jingyan Yan, Haijing Wang, Zhenyuan Cai, Jiuxiang Xie and Tongzuo Zhang in Evolutionary Bioinformatics

Click here for additional data file.

Supplemental material, sj-tif-3-evb-10.1177_1176934321996353 for Gut Microbiome Changes in Captive Plateau Zokors (Eospalax baileyi) by Daoxin Liu, Pengfei Song, Jingyan Yan, Haijing Wang, Zhenyuan Cai, Jiuxiang Xie and Tongzuo Zhang in Evolutionary Bioinformatics

Click here for additional data file.

Supplemental material, sj-tif-4-evb-10.1177_1176934321996353 for Gut Microbiome Changes in Captive Plateau Zokors (Eospalax baileyi) by Daoxin Liu, Pengfei Song, Jingyan Yan, Haijing Wang, Zhenyuan Cai, Jiuxiang Xie and Tongzuo Zhang in Evolutionary Bioinformatics