Marked Seasonal Variation in Structure and Function of Gut Microbiota in Forest and Alpine Musk Deer

Musk deer (Moschus spp.) is a globally endangered species due to excessive hunting and habitat fragmentation. Captive breeding of musk deer can efficiently relieve the hunting pressure and contribute to the conservation of the wild population and musk supply. However, its effect on the gut microbiota of musk deer is unclear. Recent studies have indicated that gut microbiota is associated with host health and its environmental adaption, influenced by many factors. Herein, high-throughput sequencing of the 16S rRNA gene was used based on 262 fecal samples from forest musk deer (M. berezovskii) (FMD) and 90 samples from alpine musk deer (M. chrysogaster) (AMD). We sought to determine whether seasonal variation can affect the structure and function of gut microbiota in musk deer. The results demonstrated that FMD and AMD had higher α-diversity of gut microbiota in the cold season than in the warm season, suggesting that season change can affect gut microbiota diversity in musk deer. Principal coordinate analysis (PCoA) also revealed significant seasonal differences in the structure and function of gut microbiota in AMD and FMD. Particularly, phyla Firmicutes and Bacteroidetes significantly dominated the 352 fecal samples from captive FMD and AMD. The relative abundance of Firmicutes and the ratio of Firmicutes to Bacteroidetes were significantly decreased in summer than in spring and substantially increased in winter than in summer. In contrast, the relative abundance of Bacteroidetes showed opposite results. Furthermore, dominant bacterial genera and main metabolic functions of gut microbiota in musk deer showed significant seasonal differences. Overall, the abundance of main gut microbiota metabolic functions in FMD was significantly higher in the cold season. WGCNA analysis indicated that OTU6606, OTU5027, OTU7522, and OTU3787 were at the core of the network and significantly related with the seasonal variation. These results indicated that the structure and function in the gut microbiota of captive musk deer vary with seasons, which is beneficial to the environmental adaptation and the digestion and metabolism of food. This study provides valuable insights into the healthy captive breeding of musk deer and future reintroduction programs to recover wild populations.


INTRODUCTION
Musk deer (Moschus spp.), the only extant genus of the family Moschidae, consists of seven species and are widely found in forests and mountains in Asia (Yang et al., 2003;Jiang et al., 2020). China has the largest musk deer population and musk source globally . Six species of genus Moschus are found in China, among which the forest musk deer (FMD, M. berezovskii) and the alpine musk deer (AMD, M. chrysogaster) are the most widely distributed and have the highest wild and captive population (Fan et al., 2018). They inhabit high-altitude coniferous forests or broad-leaved forests (Wang and Harris, 2015), alpine shrub meadow, and mountain forest grassland zone (Harris, 2016) in central and southwestern China, with some overlapping areas. However, wild musk deer populations plummeted to the brink of extinction in the 20th century due to illegal hunting, habitat fragmentation, and other human activities (Wu and Wang, 2006;Cai et al., 2020). The International Union for Conservation of Nature (IUCN) Red List (Wang and Harris, 2015;Harris, 2016) and red list of China's vertebrates (Jiang et al., 2016) have listed both the species as endangered (EN) and critically endangered (CR), respectively. The captivity breeding of the FMD and AMD began in China in the 1960s to curb the rapid decline of the musk deer population by reducing the pressure on hunting wild musk deer to some extent (Fan et al., 2019). Captive individuals can also serve as a rewilding resource for reintroduction, which is beneficial for effective conservation and population recovery of wild musk deer. However, gastrointestinal diseases occur in the captive FMD and AMD with a fatality rate of about 30%, especially in winter and autumn . Diseases are the most significant constraints on population growth, breeding scale, and musk secretion (Zhao et al., 2011;Zhang et al., 2019).
Gut microbiota has a close mutualistic symbiotic relationship with their hosts during long-term coevolution and is essential in organisms (Nicholson et al., 2012). The complex and variable micro-ecosystem of gut microorganisms were involved in metabolism, immune regulation, intestinal development promotion, pathogen defense, and other physiological activities Yoo et al., 2020). Many factors, such as host genetics, diet, season, age, and lifestyle, influence the composition and function of gut microbiota (Zhang et al., 2010;Claesson et al., 2012;O'Toole and Jeffery, 2015). For instance, changes in dietary composition rapidly alter the composition and abundance of gut microbiota (Zmora et al., 2019). Seasonal changes in food composition and availability change the structure and function of gut microbiota in many animals (Amato et al., 2015;Xue et al., 2015). Moreover, captive breeding and ex situ conservation are effective for the conservation of endangered species. Long-term captive breeding has significantly changed the dietary composition of musk deer compared with wild populations (Guo et al., 2019), resulting in changes in the composition and function of gut microbiota (Sun et al., 2020). Therefore, studying the diversity of gut microbiota of endangered species in captivity is essential for assessing the current captive conditions, understanding the appropriate capacity of changes in gut microbiota for their future, and evaluating whether they can be released into the wild. However, the relationship between the diversity, structure, and function of gut microbiota of captive musk deer of different ages and seasonal changes is unclear due to the inadequate relevant data.
In this study, combined with weighted gene co-expression network analysis (WGCNA), the 16S rRNA gene amplicon technology was used for high-throughput sequencing on the Illumina MiSeq sequencing platform to analyze the fecal microbial composition, diversity, and function in captive FMD and AMD in different seasons. This study aimed to (i) explore the composition and differences in the gut microbiota of both musk deer in different seasons; (ii) analyze the core microbiota and its metabolic functions, and their seasonal difference; and (iii) construct a weighed co-occurrence network of gut microbiota for identifying modules of co-occurring taxa and hub OTUs significantly related with the seasonal variation. Therefore, this study can provide a scientific basis for the effective management of captive musk deer.

Sample Collection
A non-invasive sampling method was used to collect 262 fresh feces of captive FMD in Aru Township, Qilian County, Qinghai Province, China (100 • 21 E, 38 • 7 N), during the early spring (mid-March, here referred to as T1, N = 49), late spring (late May, here referred to as T2, N = 57), summer (mid-July, here referred to as T3, N = 56), autumn (mid-November, here referred to as T4, N = 50), and winter (late December, here referred to as T5, N = 50). The FMD breeding center was located northeast of the Qinghai-Tibet Plateau with an altitude, annual average temperature, and precipitation of about 3,002 m, −0.1 • C, and 403 mm, respectively, with high daily range and radiation intensity. The maximum and minimum temperatures appeared in July and January, respectively. Furthermore, 90 fecal samples of captive AMD were collected during the late spring (late May, here referred to as T2, N = 55) and winter (late December, here referred to as T5, N = 35) in Xinglong Mountain, Gansu Province (104 • 4 E, 35 • 49 N). The AMD breeding center was located northeast of the Qinghai-Tibet Plateau with an altitude, annual average temperature, and precipitation of 2,171 m, 5.4 • C, and 406 mm, respectively.
The FMD and AMD breeding centers were cleaned the night before sampling, and each individual was kept in a separate enclosure so that the fresh feces of each individual could be collected the following morning. Disposable PE gloves were used to collect fresh feces and put them into sterile self-sealed bags to prevent contamination of the feces surface. All samples were temporarily stored in the −20 • C vehicle-mounted refrigerator and later stored in the −80 • C ultra-low temperature refrigerator in the laboratory for DNA extraction.
DNA Extraction and 16S rRNA Gene Sequencing E.Z.N.A. R soil DNA Kit (Omega Biotek, Norcross, GA, United States) was used to extract total bacterial DNA from FMD and AMD feces. The extraction quality of DNA was assessed using 1% agarose gel electrophoresis. NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, United States) was used to determine the concentration and purity of DNA. The V4-V5 variable region of the 16S rRNA gene was amplified via PCR using primers; 515F (5 -GTGCCAGCMGCCGCGG-3 ) and 907R (5 -CCGTCAATTCMTTTRAGTTT-3 ). Each sample had three PCR replicates.
The PCR was run in a reaction volume of 20 µl, containing 4 µl TransStart FastPfu buffer (5×), 2 µl dNTPs (2.5 mM), 0.8 µl each of forward and reverse primers (5 µM), 0.4 µl TransStart FastPfu DNA Polymerase, and 10 ng sample DNA and topped up with ddH 2 O. Throughout the PCR amplification, ultrapure water was used instead of a sample solution as a negative control to eliminate the possibility of false-positive PCR results. The PCR amplification procedure was as follows: 95 • C for 3 min (initial denaturation), followed by 27 cycles at 95 • C for 30 s (denaturing), 55 • C for 30 s (annealing), and 72 • C for 45 s (extension), and a final single extension at 72 • C for 10 min.
The same sample PCR products were mixed, then 2% agarose gel was used to recover the PCR products. AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, United States) was used to purify the recovered products. The recovered products were measured using 2% agarose gel electrophoresis and quantified using Quantus TM Fluorometer (Promega, United States). NEXTFLEX Rapid DNA-Seq Kit (Bioo Scientific, Austin, TX, United States) was used to build the library according to the manufacturer's protocols. Then, the purified PCR-amplified fragments were pooled in equal concentrations and sequenced on Illumina MiSeq PE300 platform (San Diego, CA, United States).

Determination of OTU and Taxonomy Assignments
The raw data were pre-processed to remove the known adaptor, specific primers and low-quality ends with Trimmomatic (version 0.39, PE-phred33 ILLUMINACLIP:2:30:10 TRAILING:20 MINLEN:50 SLIDINGWINDOW:50:20) (Bolger et al., 2014). We set a 50-bp slide window and trimmed off those sequences with average base quality <20. These lowquality reads include reads with >10 nucleotides aligned to the adapter sequences, those with unidentified nucleotide (N) sequences, and those read lengths below 50 bp. The generated forward and reverse unpaired sequences were than merged together using FLASH with a minimum overlap of 10 bp and maximum mismatch of 0.2 (version 1.2.7, -m 10 -x 0.2) (Magoč and Salzberg, 2011).
Chimera reads were removed, and operational taxonomic units (OTUs) were clustered with 97% nucleotide sequence similarity using UPARSE (version 7.1) 1 (Costello et al., 2009). The highest-frequency sequence in each OTU was selected as the representative sequence for further annotation. The reference sequence annotation and curation pipeline (RESCRIPt) was used to prepare a compatible Silva 138/16s bacteria database, based on SILVA SSURef of the curated NR99 (version 138) database 2 followed the protocol suggested by author 3 . The Ribosomal Database Project (RDP) classifier (version 2.11) 4 was used to classify the representative OTU sequence against the Silva 138/16s bacteria database at a confidence threshold of 0.8 . The taxonomy-based filtering was applied to remove all features that contain either mitochondria, chloroplast, or archaea, and the results were aligned to generate the OTU table.

Bioinformatic Analysis
Operational taxonomic units that were present in any two samples that are less than five sequences as well as a total abundance (summed across all samples) of less than 10 also be filtered. The OTU table was than rarefied to the lowest number of reads across samples (for FMD of 49,615 and for AMD of 50,680) before downstream analysis. The corresponding abundance information of each OTU annotation result in each sample was counted, and the sample sequences were defined as OTU based on the minimum number of sample sequences. Community bar charts were used to plot the relative abundance of bacteria in each fecal sample in musk deer at phylum and genus levels with R software (version 3.3.1, packages "stats"). Venn charts were used to analyze the core and unique bacterial phyla and genera in different seasons with the R software. A cluster heat map was used to compare the composition in different seasons with the R software (packages "pheatmap") (Perry, 2016). Alpha diversity was used to reflect the diversity of intestinal microbial composition. Sobs index (the observed richness) and Shannon index were used to analyze the diversity of gut microbial composition at the OTU level using Qiime software 5 (Caporaso et al., 2010). The Wilcoxon rank-sum test was used to analyze the seasonal significance of the two Alpha diversity indexes with the R software (packages "stats").
The bacterial diversity among different microbial communities was compared and analyzed to explore the community composition among different seasonal groups. Principal coordinate analysis (PCoA) was used to analyze the beta diversity among different groups with the R software (packages "vegan"). Bray-Curtis was used to calculate the distance between samples at the OTU level. Analysis of similarities (ANOSIM), a non-parametric statistical test, was used for the intergroup difference test with the R software (packages "vegan, " anosim function) (Oksanen et al., 2019). The metabolic functions of bacterial communities were predicted using PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) software based on the KEGG (Kyoto Encyclopedia of Genes and Genomes) and EggNOG (Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups) databases with the OTU species annotation and abundance information (Langille et al., 2013). Moreover, the Wilcoxon rank-sum test was used to analyze the seasonal differences of metabolism-related functions and dominant bacteria in all groups.
Weighted gene co-expression network analysis is a system biology method originally conceived for describing correlation patterns among genes across microarray data, which is widely used to identify critical hub genes of biological processes by constructing gene co-expression networks (Saris et al., 2009). Currently, it has been also applied in gastric microbiome networks and microbial modules (Park et al., 2019). We used WGCNA analysis to identify modules of co-occurred taxa and relate these modules to traits of musk deer (season, age, and gender). R package WGCNA (version 1.70-3) 6 was used to perform unsigned WGCNA analysis (Langfelder and Horvath, 2008). The weighted correlation network analysis was performed to cluster OTUs (van Dam et al., 2018).
The soft thresholding power of 8 was chosen based on the criterion of approximate scale-free topology as well as the mean connectivity lower than 100 (Morandin et al., 2016). To divide highly co-occurred OTUs into several module members, hierarchical clustering using the dynamic tree cut method with DeepSplit of 3 was performed to create a hierarchical clustering tree of OTUs as a dendrogram (Do et al., 2017). Module eigengenes were calculated using the moduleEigengenes() function in the WGCNA R package to demonstrate the correlation of the module eigenvalue and OTU abundance profile . Then, we identified potential modules and OTUs associated with season, age, and gender. The network was generated to visualize the correlations among OTUs in the module associated with season . The connection weight values (range from 0 to 1) were calculated using the intramodularConnectivity() function in the WGCNA R package to indicate the strength of co-regulation between taxa. Additionally, Cytoscape software (version 3.8.2) 7 was used for network visualization, with Cytohubba plugin analyzing and extracting the hub OTUs (Kavarthapu et al., 2020).

Assessment of Sequence Data
A total of 36,862,468 (140,696 reads/sample) and 12,912,306 (143,470 reads/sample) high-quality clean reads were obtained in FMD and the AMD samples, respectively. The rarefaction curves of sobs and Shannon indexes at the OTU levels became gradually placid as the sequencing depth increased (Supplementary Figure 1). The results demonstrated that each fecal sample had sufficient OTUs to reflect the maximum level of bacterial diversity, which indicated a sufficient sequencing depth.
A total of 3,548 and 2,259 OTUs were identified in FMD and AMD, respectively. At a 97% sequence identity threshold, the OTUs of the FMD were classified into 20 phyla, 35 classes, 93 orders, 171 families, and 404 genera, while those of the AMD were classified into 17 phyla, 27 classes, 66 orders, 125 families, and 300 genera.

Seasonal Variation in α Diversity
The Good's Coverage index of captive FMD and AMD in different seasons was higher than 99%. Both sobs and Shannon indexes were used to reflect the diversity of gut microbiota of musk deer in different seasons.
For captive juvenile FMD, the α diversity of gut microbiota was significantly lower in early spring than in other periods (p < 0.05), and it was higher in winter than in other periods (Figures 1G,H). For adult FMD, the α diversity was significantly higher in winter than in other periods (p < 0.05) and lower in summer than in other periods. Overall, the α diversity of gut microbiota in captive FMD was higher in winter and autumn than in spring and summer.
For captive adult AMD, the α diversity of gut microbiota was significantly higher in winter than in late spring (Figures 1G,H) FIGURE 1 | Gut microbial composition of the FMD and AMD. Histogram of relative abundance of individual bacteria phyla of the FMD (A) and the AMD (G). Relative abundance of dominant phyla in the FMD (A) and the AMD (B). Heat map analysis based on identifiable bacterial genera with the relative abundance of top30 for the FMD (C) and the AMD (D). The red, blue, light green, black, and dark green letters represent phylum Firmicutes, Bacteroidetes, Planctomycetes, Proteobacteria, and Spirochaetes, respectively. Analysis of core and unique bacteria of the FMD (E) and AMD (F) at phylum and genus levels via Venn plots. Seasonal variation of α diversity in the gut microbiota of the musk deer based on Sobs index (G) and Shannon index (H). * p < 0.05 (Wilcoxon rank-sum test), * * p < 0.01, and * * * p < 0.001. ns, not significant.
(p < 0.05). Besides, the α diversity of gut microbiota was higher in winter than, but not significant, for both juvenile and older AMD.

Seasonal Variation in β Diversity
The Bray-Curtis distance algorithm was used to calculate the distance between fecal samples of captive musk deer in different seasons. ANOSIM analysis was used to test whether the intergroup differences were significantly greater than the intra-group differences in different seasons. The PCoA analysis showed that the R values were all greater than 0 (p = 0.001).
Besides, there were significant seasonal differences in the gut microbial composition of the FMD in different seasons, and the intergroup differences were significantly greater than the intragroup differences (Figures 2A-C). The significant differences between the summer and winter groups (R = 0.3078, p = 0.0010) and the autumn and winter groups (R = 0.3254, p = 0.001) were relatively high. The gut microbial composition of captive adult AMD was significantly different between the winter and late spring groups (R = 0.1332, p = 0.016). Moreover, the intergroup difference was significantly greater than the intragroup difference ( Figure 2E). However, there was no significant seasonal difference in the composition of gut microbiota in juvenile (R = 0.0392, P = 0.272) or older (R = −0.0067, P = 0.512) AMD (Figures 2D,F).

Analysis of Seasonal Difference of Dominant Bacteria
The Wilcoxon rank-sum test was used to analyze the significant seasonal differences. For both juvenile and adult FMD, the  Figures 2, 3). PCoA analysis between late spring and winter groups in juveniles (D), adults (E), and older (F) AMD.
relative abundance of Firmicutes was significantly higher in winter than in other periods (p < 0.05) (Figure 3A), while that of the Bacteroidetes showed the opposite result (p < 0.05) ( Figure 3B). The relative abundance of Firmicutes in summer was lower than that in other periods, while that of Bacteroidetes showed the opposite result. Overall, the relative abundance of Firmicutes in the FMD was higher in autumn, winter, and early spring than in late spring and summer. In contrast, that of Bacteroidetes showed the opposite result (Figures 3A,B). Moreover, for captive AMD of different ages, the relative abundance of Firmicutes was higher in winter than in late spring, while that of Bacteroidetes showed the opposite result. Particularly, there were significant seasonal differences among adult AMD (p < 0.05). The analysis of the ratio of Firmicutes to Bacteroidetes (F/B ratio) showed that the ratio was significantly higher in winter than in other periods for both juvenile and adult FMD (p < 0.05) (Figures 3C,D). Overall, the F/B ratio significantly decreased in summer than in spring and increased in winter than in summer. Similarly, for captive AMD of different ages, the F/B ratio was higher in winter than in late spring, and there were significant seasonal differences among adult AMD (p < 0.05) (Figure 3E).
The seasonal difference analysis indicated that nine dominant identifiable bacteria genera in juvenile and adult FMD had significant seasonal differences. The relative abundances of the genera Christensenellaceae R7 group and Ruminococcus were significantly higher in winter and early spring than in other periods. The relative abundances of the genera UCG-005 and Monoglobus were significantly higher in summer than in other periods. The relative abundance of genus Alistipes was significantly higher in late spring and summer than in other periods (Figures 3F,G). Moreover, among the nine dominant bacteria genera in AMD, genera Bacteroides, UCG-005, Alistipes, and Anaerostipes showed significant seasonal differences. The relative abundances of genera UCG-005 and Alistipes were significantly higher in winter than in late spring, while those of Bacteroides and Anaerostipes showed opposite results ( Figure 3H).

Functional Prediction Analysis
The functional prediction analysis in the KEGG database showed that gut microbiota in captive FMD and AMD were mainly involved in carbohydrate metabolism, amino acid metabolism, energy metabolism, metabolism of cofactors and vitamins, metabolism of cofactors and vitamins, translation, and replication and repair. The EggNOG database showed that 17 functions had high relative abundance.
The KEGG database showed that the four major metabolismrelated functions of captive FMD had seasonal differences. Overall, those functions were significantly lower in summer than in other periods (p < 0.05) and were relatively higher in spring than in other periods ( Figure 4A). The main metabolism-related function of captive AMD was higher in late spring than in winter. Also, the energy metabolic function had significant seasonal differences (p < 0.05) (Figure 4C).
The EggNOG database showed that six major metabolismrelated functions also had seasonal differences. Overall, the functions were relatively higher in spring than in other periods. All the major metabolic functions of gut microflora in captive FMD were significantly lower in summer than in other periods except for lipid transfer and metabolic function (p < 0.05) (Figure 4B).
Similarly, all the major metabolic functions in captive AMD were higher in late spring than in winter except for amino acid transfer and metabolic function ( Figure 4D).

The Hub OTUs Fluctuated With the Season Using WGCNA Analysis
Baseline characteristics of participants in weighted correlation network analysis. A total of 262 FMD and 3493 OTUs were FIGURE 4 | Gut microbial function prediction and seasonal difference analysis based on KEGG and EggNOG databases. Seasonal difference analysis of the main function of gut microbiota of the FMD based on the KEGG database (A) and EggNOG database (B). Seasonal difference analysis of the main function of gut microbiota of the AMD based on the EggNOG database (C) and KEGG database (D). * p < 0.05 (Wilcoxon rank-sum test), * * p < 0.01 and * * * p < 0.001. ns, not significant.
included in the study, with the mean age and the proportion of males of about 2.45 and 75.19%, respectively. The individual in early spring (T1), late spring (T2), summer (T3), autumn (T4), and winter (T5) periods accounted for 18. 70, 21.76, 21.37, 19.08, and 19.08%, respectively. The sample dendrogram and trait heat map is shown in Supplementary Figure 4.
The network was constructed using one-step network construction. The networkType was set to signed and softthreshold power 8 to define the adjacency matrix based on the criterion of approximate scale-free topology (Supplementary Figure 5), with minimum module size 30, the module detection sensitivity DeepSplit set to 3, and cut height for merging of modules of 0.25 which means that the modules whose eigengenes are correlated above 0.75 will be merged (Supplementary Figure 6).
A total of 3,493 OTUs were parsed into 13 different color modules. Among these 13 modules, the gray module indicated unassigned bacterial taxa. In the dendrogram, each leaf, represented as a short vertical line, corresponded to a bacterial taxon. Densely interconnected branches of the dendrogram group represented highly co-occurring bacterial taxa (Supplementary Figure 7).
The correlation between module eigenvalue and trait was calculated. The module-trait relationship heat map demonstrated the correlation coefficient between module eigenvalues and traits (−1 to 1). The black module was significantly correlated with five different periods simultaneously (p < 0.01) and was not correlated with other traits (age and gender) ( Figure 5A). Therefore, the black module was selected for subsequent analysis.
Furthermore, the OTU in the top 10% of global importance [Ktotal (the global network connectivity) value accounted for the top 10%] in WGCNA was extracted, and the network diagram was constructed in Cytoscape (715 nodes and 9,447 edges). The results showed that the black module was relatively independent and had low correlation with other modules (Figure 5B). The black module had the highest correlation with seasonal change, and 14 OTUs in the global hub OTUs were located in the black module, which was significantly related to seasonal changes ( Figure 5C). OTUs with a connectivity weight value greater than 0.02 in black modules were extracted to construct the network diagram of black modules (Figure 5D). The thickness of gray lines represented the co-occurred strength between OTUs. The size of an OTU node represented the Kwithin (the intramodular connectivity) value and the connectivity between the node and the black module. CytoHubba plugin was further used to extract and identify the hub OTUs including OTU6606, OTU5027, OTU7522, and OTU3787 ( Figure 5E) as top 4, which belong to the order Oscillospirales in phylum Firmicutes, and the first three hub OTUs belong to genus UCG-005 of family Oscillospiraceae.

DISCUSSION
The normal steady-state gut microbiota is closely related to the health of the host and coevolves in a mutualistic relationship with the host through complex interactions. Its unique community structure and metabolites are essential for regulating host metabolism, growth, and development, resistance to pathogens, immune regulation, adaptation, evolution, etc. (Bäckhed et al., 2005;Engel et al., 2012;Ezenwa et al., 2012;Zhang et al., 2016). This study systematically and comprehensively analyzed the seasonal differences of gut microbiota in captive FMD and AMD at different ages. The results indicated that the FMD and AMD shared seven dominant bacteria genera, including the genera Bacteroides, Christensenellaceae R7 group, UCG 005, Rikenellaceae RC9 gut group, Alistipes, Ruminococcus, and Prevotellaceae UCG 004 and genera Christensenellaceae R7 group, Ruminococcus, and Prevotellaceae UCG 004. Current studies have shown that genus Christensenellaceae R7 group is widely found in the intestinal tract and mucosa of the host, which is essential for good health and is involved in amino acid and lipid metabolisms (Waters and Ley, 2019). Genus Ruminococcus is involved in the degradation of cellulose and hemicellulose in the rumen of ruminants (Matulova et al., 2008;La Reau et al., 2016). It produces several cellulases and hemicellulases that convert dietary fiber into various nutrients (La Reau and Suen, 2018). It is also involved in food digestion and carbohydrate metabolism in ruminants. Genus Prevotellaceae UCG 004 is involved in the degradation of polysaccharides and the production of short-chain fatty acids (SCFAs) (Heinritz et al., 2016). Family Christensenellaceae (genus Christensenellaceae R7 group belongs to this family), Ruminococcaceae (Ruminococcus belongs to this family), and genus Alistipes are associated with immune regulation and healthy homeostasis and are regarded as potentially beneficial bacteria (Kong et al., 2016;Shang et al., 2016;Wang et al., 2018). Among them, the genera Christensenellaceae R7 group and Alistipes are classified as potential biomarkers for intestinal diseases (Crohn's disease, colorectal cancer, Clostridium difficile infection, etc.) (Mancabelli et al., 2017). Moreover, the genera Bacteroides, Alistipes, and Rikenellaceae RC9 gut group belong to Bacteroidetes. Genus Bacteroides can improve the metabolism of the host, mainly by participating in the metabolism of bile acid, protein, and fat, and regulating carbohydrate metabolism. Genus Alistipes is involved in the metabolism of SCFAs. Bacteroides and Alistipes belong to bile-tolerant microorganisms (David et al., 2014). High-fat animal feed can increase the relative abundance of genera Bacteroides and Alistipes (Wan et al., 2019), thus improving lipid metabolism by regulating acetic acid production (Yin et al., 2018). Genus Rikenellaceae RC9 gut group also promotes lipid metabolism (Zhou et al., 2018).
Furthermore, phyla Firmicutes and Bacteroidetes significantly dominated the 352 fecal samples from captive musk deer in different seasons (relative abundance of more than 90%), consistent with other studies Zhao et al., 2019;Fountain-Jones et al., 2020). Firmicutes and Bacteroidetes are the dominant core bacteria in the rumen of ruminants, with the highest relative abundance. They are involved in essential processes, such as food digestion, nutrient regulation, and absorption, energy metabolism, and defense against invasion of foreign pathogens in the gastrointestinal tract of hosts (Jewell et al., 2015;Guan et al., 2017;Wang et al., 2017). Firmicutes promote fiber degradation in food and convert cellulose into volatile fatty acids, enhancing food digestion and growth and development. Herein, among the top 30 bacterial genera of gut microbiota, 23 and 19 were Firmicutes in the FMD and AMD, respectively. Besides, most were associated with carbohydrate metabolism and cellulose digestion and absorption. Bacteroidetes mainly promote the digestion and absorption of proteins and carbohydrates in the food, thus promoting the development of the gastrointestinal immune system. Among the top 30 identifiable bacterial genera, five and eight were Bacteroidetes in the FMD and AMD, respectively, of which most were involved in lipid metabolism. The nutrient level of animal feed directly affects the abundance of Firmicutes and Bacteroidetes. Herein, Firmicutes and Bacteroides were the core dominant microflora in the FMD and AMD at different seasons and ages. In addition, the relative abundances of Firmicutes and Bacteroidetes and their ratios had significant seasonal differences. The relative abundances of Firmicutes and the F/B ratio were higher in the cold season than in the warm season. This is beneficial for captive FMD and AMD to adapt to the cold season, thus promoting the decomposition of cellulose and hemicellulose in feed, carbohydrate metabolism, and nutrient digestion and absorption.
The FMD and AMD are typical ruminants. Wild individuals feed on various high-fiber leaves, while captive individuals feed on concentrated animal feed and rouge feed. The former mainly includes carrot (Daucus carota), potato (Solanum tuberosum), maize (Zea mays), soybean (Glycine max), etc., and the latter includes fresh leaves (warm season) or dry leaves (cold season). Core and dominant bacteria genera are essential in food digestion, nutrient absorption, and energy metabolism of captive FMD and AMD. Herein, the relative abundances of the genera Christensenellaceae R7 group, Ruminococcus, and Prevotellaceae UCG 004 (Firmicutes) were higher in the cold season than in the warm season. In contrast, the relative abundance of genus Bacteroides (Bacteroidetes) showed the opposite results, indicating that the seasonal difference of dominant bacteria genera is closely associated with animal feed composition in different periods.
Principal coordinate analysis showed that both adult and juvenile musk deer had significant seasonal differences, indicating that the gut microbial structure and composition of FMD and AMD are significantly different at different seasons, consistent with other studies (Amato et al., 2015;Maurice et al., 2015;Trosvik et al., 2018). The seasonal difference is closely related to food resources, dietary structure, nutrient utilization, and feeding pattern (Xue et al., 2015;Wu et al., 2017;Kartzinel et al., 2019). The Wilcoxon ranksum test showed that the α diversity of the gut microbiota had significant seasonal differences in captive FMD and AMD. Overall, the α diversity of the gut microbiota in musk deer was higher in winter and autumn than in spring and summer. Previous studies have shown that higher α diversity leads to a more complex and stable intestinal microbiota composition, enhancing resistance to external interference, and adaptability, which is beneficial to the host health (Stoffel et al., 2020). The decrease or loss of α diversity is associated with various diseases (Rogers et al., 2016). Therefore, the reported higher α-diversity of the gut microbiota in captive FMD and AMD in autumn and winter can improve the resistance to adverse environmental factors, reduce the influence of adverse environmental factors, and promote the intake of fiber-rich food and nutrient absorption and utilization in the cold season.
The function prediction analysis based on KEGG and EggNOG databases showed that the gut microorganisms in captive FMD and AMD mainly involve carbohydrate metabolism, amino acid metabolism, energy metabolism, and cofactor and vitamin metabolism. These metabolic functions also had seasonal differences and were significantly correlated with core bacteria phyla and genera. The seasonal changes of core bacteria were also closely associated with the seasonal function changes. The effect of season variation on gut microbiota of musk deer was higher than that of age and gender based on the correlation between module eigenvalue and traits of musk deer. WGCNA analysis indicated that the black module had low correlation with other modules and has the highest correlation with seasons. OTU6606, OTU5027, OTU7522, and OTU3787 were significantly fluctuated with the season, which all belong to phyla Firmicutes. The four hub OTUs were at the core of the black module and can greatly affect the network structure of the co-occurrence bacterial taxa network in the black module, which can be used as an important indicator to evaluate the intestinal health of captive forest musk deer.

CONCLUSION
The study systematically and comprehensively analyzed the seasonal differences of gut microbiota structure and function through 16S rRNA gene sequencing of 352 fecal samples from the captive FMD and AMD at different ages. Comparison analysis identified significant seasonal variations of α-diversity, core microbiota, and metabolic functions. The α-diversity, phylum Firmicutes, and F/B ratio in the gut microbiota of both musk deer were higher in the cold season than in the warm season. The major metabolism-related functions were also significantly higher in the cold season than in the warm season in FMD. The four identified hub OTUs can be used as an important indicator to evaluate the intestinal health of captive forest musk deer. Therefore, this study provides insights into the captive breeding environment and future reintroduction programs. Besides, functional annotation analysis of gut microbiota and the seasonal changes of metabolic pathways combined with metagenome and metabolomic analyses can help to further explore the role of gut microbiota in the health, environmental adaptation, and metabolism of the FMD and AMD in the future.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the NCBI GenBank, accession number PRJNA725631 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA725631?review er=e2rcdic27nir1a7qhp8unlk9s6).

ETHICS STATEMENT
A total of 262 fresh feces of captive forest musk deer were collected with the permission of the authorities of the forest musk deer breeding center in Qilian county, China. Also, a total of 90 fecal samples of captive alpine musk deer were collected with the permission of the authorities of the alpine musk deer breeding center in Yuzhong county, China. All procedures were approved by the Ethical Committee for Experimental Animal Welfare of the Northwest Institute of Plateau Biology.