Quality filtering and validation of the dataset
There was an average of 909,993 reads (average length = 428.67 bp; Table S1) obtained from each sample after quality filtering. There were 1,693 OTUs detected, and each sample had between 396 and 865 OTUs (Table S2). Individually rarefaction curves showed that the sequencing depth was adequate for the study (Fig S1). Good’s coverage showed that 28 samples had between 99.28 to 99.64% coverage (Table S2). The results showed that the number of OTUs increased rapidly until it gradually reached a plateau. It shows that the amount of sequencing data is large enough, and the increase of sequencing data can not find more OTU, which can reflect the information of most microbial species in the sample.
Core bacterial communities and composition difference in the two species
There were 21 phyla, 34 classes, 82 orders, 142 families, and 310 genera identified in the two musk deer species. There were 658 OTUs that were present in all groups. There were 100, 52, 128, and 10 OTUs that appeared in each sample of the FS, FW, SS, and SW groups, respectively (Fig 1A; FMD in summer and winter: FS and FW; SMD in summer and winter: SS and SW ).
The classification of sequences from the samples resulted in five different phyla was identified, including Firmicutes, Bacteroidetes, Spirochaetes, Verrucomicrobia, and Proteobacteria (Fig 1B). We considered these taxa that were shared by all individuals in a group to be the core bacterial community for that group. Firmicutes (FMD: 53±8.31 to 64.95±11.73, SMD: 86.2±5.713 to 90.21±3.265) and Bacteroidetes (FMD: 38.06±7.431 to 29.03±9.044, SMD:12.75±6.008 to 8.876±2.31) were core and dominant community members for both musk deer species. The relative abundances of Firmicutes were significantly higher in SMD than in FMD, and this corresponded to decreases in the relative abundances of Bacteroidetes (Fig 1B；Table 2). Spirochaetes and Verrucomicrobia were only found in forest musk deer (Fig 1B). Also, the difference in gut microbiota composition between the two musk deer can be clearly observed according to the Heatmap . Our results indicated the existence of two core taxa at the class level (Clostidiales and Bacteroidales) (Fig S2). At the order level, core gut bacterial included Bacteroidaceae, Lachnospiraceae, Oscillospirales, and Christensenellaceae (Fig S3), in addition to this 11 core taxa at the family level(Fig S4). These results indicated that the composition of the microbial community may be related to different physiological characteristics in the gut of musk deer.
The alpha diversity of forest musk deer was higher than that of Siberian musk deer based on the Shannon and Chao indices (Shannon index: FMD=4.682±0.428, SMD=4.54±0.56, p=0.431; Chao index: FMD=826.46±133.04, SMD=764.18±144.89, p=0.347; Table S3;Fig 2A-2B ), it is concluded that the community abundance and community diversity of the forest musk deer are higher than the siberian musk deer, but there is no significant difference. There were differences in the beta diversity of the SMD and FMD communities PCoA analysis (Fig 3). ANOSIM analysis reveals intergroup differences in bacterial communities between two musk deer species (R = 0.6655, p = 0.001). The results of ANOSIM analysis and PCoA analysis mutually verified the differences of bacterial communities between the two musk deer species.
Seasonal differences in bacterial communities
The alpha diversity of the bacterial community in the summer was higher than in the winter for the two musk deer species(Shannon index; FS:4.913±0.249, SS: 4.857±0.517, FW: 4.451±0.454, SW: 4.222±0.443). There was a significant difference in Shannon diversity index between the two seasons for FMD (p = 0.026,Q=0.031), but not for SMD ( p = 0.194, Q=0.233; Fig 2C; Table S3). The Chao index have significant different in SMD (Chao:p = 0.03, Q=0.061; Fig 2D). At the phyla level, there were 9 taxa (Firmicutes, Bacteroidota, Verrucomicrobiota, Actinobacteriota, Cyanobacteria, Patescibacteria, Elusimicrobiota, Myxococcota, and Deferribacterota )whose relative abundance was significantly different between seasons in FMD; at the same time, and only one taxon (unclassified_k__norank_d__Bacteria) whose relative abundance was significantly different between seasons in SMD (Table 2; Fig S5 ). At the genus level, there were 14 taxa (UCG-005, norank_f__norank_o__RF39, Phascolarctobacterium, NK4A214_group, Papillibacter, UCG-002, Eubacterium_siraeum_group, Arthrobacter, unclassified_k__norank_d__Bacteria, norank_f__norank_o__Veillonellales-Selenomonadales, Akkermansia, Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium, Brevundimonas, and unclassified_f__Comamonadacea )whose relative abundance was significantly different between seasons in SMD(Fig 4A). It is worth mentioning that the Firmicutes/Bacteroidetes ratio changed with the season in two musk deer.
We are much interested in musk deer gut microbiota is to determine whether specific microbiota could be the change with different season factors. To identify specific microbial candidates present in two specise seasons, LEFse (LDA score ≥ 4) analysis were used to determine specific biomarker belonged to different taxonomy (Fig 4B) from phylum to genus. The result revealed that twelve families (Bacteroidia, Bacteroidales, Bacteroidetes, Prevotellaceae, Prevotellaceae_UCG_004, Bacteroides, Bacteroidaceae, Verrucomicrobia, Verrucomicrobiae, Verrucomicrobiales, Akkermansia, and Akkermansiaceae) were significantly enriched in the FS group, three families (Lachnospirales,Lachnospiraceae, and unclassified Lachnospiraceae) were significantly enriched in the FW group, five families (UCG_010, Monoglobus, Monoglobales and Monoglobaceae) were significantly enriched in the SS group, and five families (Firmicutes, Clostridia, Christensenellaceae_R_7_group, Christensenellaceae, and Christensenellales) were significantly enriched in the SW group (Fig 4B).
Predicted functions of the microbial communities
There were significant differences in the functional capacity of SMD and FMD microbial communities. Specifically, three KEGG pathways (Environmental Information Processing, Cellular Processes, and None) were enriched in SMD samples, and five KEGG pathways (Human Diseases, Organismal Systems, Metabolism, Genetic Information Processing, and Unclassified) were enriched in FMD samples at Level 1 (Fig 5A). Also, two groups of samples were significant differences in thirty-three KEGG pathways (Level 2; p < 0.005 ). Among them, seven KEGG pathways (Transcription, Membrane Transport, Cell Motility, Signal Transduction, Cellular Processes and Signaling, Environmental Adaptation, and Carbohydrate Metabolism) were significantly enriched in SMD samples. Twenty-six KEGG pathways (Translation, Replication and Repair, Nucleotide Metabolism, "Folding, Sorting and Degradation", Enzyme Families, and Metabolic Diseases, et al) were significantly enriched in FMD samples (Fig 5B).At level 3, a functional comparison of 143 extremely significant differences (p < 0.001) was performed. Among them, the SMD group enriches thirty-five functions, and the FMD group enriches other functions. The results also reflect differences in transcription and metabolism between the two groups of animals (Fig S6).
There was a significant difference in predicted microbiota function in SMD between summer and winter. The relative abundance of genes (Level 1) involved in Metabolism, Organismal Systems, and Human Diseases was significantly different in two seasons (Fig 5C). The relative abundance of gene families(Level 2) involved in Xenobiotics Biodegradation and Metabolism, Endocrine System, and Metabolism of Cofactors and Vitamins was significantly higher in the SS samples than in the SW samples (Fig 5C). Besides, the Influenza A, Toxoplasmosis, Colorectal cancer, Viral myocarditis, and Smell cell lung cancer functions (Level 3) were significantly enriched in the SS group (Fig 5C).