Taxonomic analysis
The unigene sequences were aligned using DIAMOND software. To visualize the species annotation results and effectively present the relative abundance of species at various taxonomic levels in each sample, we employed Krona. The results indicated that in all 15 samples MC2101, MC2102, MC2103, MC32013, MC3201, MC3202, MM1401, MM1402, MM1403, MM2201, MM2202, MM2203, MM3101, MM3102 and MM3103, bacteria were the most abundant microorganisms at 72-82%. They were followed by unknown microorganisms, archaea, viruses and eukaryota, regardless of the yak’s region.
Relative abundance of most abundant species at different taxonomic level
Bacillota, Bacteroidota, Actinomycetota and Euryarchaeota were more abundant phyla in all 5 groups G1, G2, G3, G4 and G5 (Fig 1A). Bacillota were more abundant in group G1, G2, G3 followed by Bacteroidota and Actinomycetota and Euryarchaeota, whereas in group G4 Actinomycetota were more abundant followed by Bacillota, Euryarchaeota and Bacteroidota and in group G5 Bacillota were more abundant followed by Bacteroidota, Euryarchaeota, Actinomycetota. At genus level, in group G4 Arthrobacter were more abundant followed by Methanobrevibacter, Peribacillus, Psychrobacillus, clostridium and Bacteroides. Whereas in group G5 Methanobrevibacter were more abundant followed by clostridium, Bacteroides, Arthrobacter and Romboutsia and in group G3 Arthrobacter were more abundant followed by Bacteroides, Alistipes, Clostridium, Methanobrevibacter, Akkermansia, Pseudoflavonifractor, in G1 Romboutsia were more abundant followed by clostridium, Bacteroidota, Alistipes, in group G2 Bacteroides were more abundant followed by Clostridium, Methanobrevibacter, Alistipes (Fig 1B). The heat map generated based on the abundance information of the top 35 genera revealed significant differences in the composition of the gut microbiota among different region groups of yaks. Each group from different was distinctly clustered from phylum to species level, suggesting region specific gut microbiota signatures (Fig 2).
The PCA results indicate that microbiota composition between different region group is different from phylum to species level (Fig 3). Further the results of principal co-ordinates analysis based on the Bray-Curtis distance also shows that microbiota community structure and composition was different between groups and was more similar within the groups from phylum to species level (Fig 4). These findings were verified by ANOSIM analysis which indicates that microbiota composition difference at different taxonomic level between groups of different Yak region was significantly greater than within group (Fig 5).
Common functional database analysis
The metagenomic sequencing profiling data of samples were annotated to the KEGG database to evaluate the functional diversity. Fig 6 and 7 show the differences in the potential functional abundance of the bacterial community between the groups. The results of the pathway annotation are provided in Fig 8 and 9. Moreover, the CAZy database analysis displayed highest abundances of glycoside hydrolases (124491) followed by glycosyl transferases (70228), carbohydrate-binding modules (31399), carbohydrate esterases (11836) and polysaccharide lyase (2366), (Fig 10 and 11).
Antibiotics resistance gene analysis
Based on the abundance table resistance genes, the content and percentage (Fig 12) of antibiotic resistance ontology (ARO) in each sample and screen out the top 20 AROs with the highest abundance was determoneed. It has seen that there was a difference in the relative abundance and percentage of ARO among different yak age groups. Among selected top 20 AROs, the relative abundance of vancomycin resistant genes was found to be highest in all groups followed by macrolide, tetracycline, oxazolidinone, bacitracin, aminoglycoside. However, vanW_gene_in_vanI_cluster was found to be highest in G1 and lowest in G4 which contains highest relative abundance of vanY_gene_in_vanB_cluster. Further clustered heat map results showed that different age yak groups samples were clustered individually and ANOSIM analysis based on the Bray-Curtis distance also showed that difference between groups was significantly higher than within groups of different age.
The observed high prevalence of unknown microorganisms in samples underlines the necessity of enhancing current microbial databases to comprehensively cover global microbial diversity
(Almeida et al., 2019). The results also revealed a varied representation of different phyla across different sample groups, with Bacillota, Bacteroidota, Actinomycetota and Euryarchaeota being the most abundant. Interestingly, these phyla’s representation did not remain consistent across all groups. Such variations in microbial compositions across different sample groups have been reported in other studies and could be attributed to various factors like diet, environment, host genetics and age
(Zhao et al., 2015). At a deeper taxonomic level, the genus Arthrobacter was found to be more abundant in groups G4 and G3, a finding that has been mirrored in previous works
(Mhuireach et al., 2016). Arthrobacter spp. are known for their remarkable metabolic versatility and resilience to environmental stressors, which might explain their abundant representation in male of older ages
(Mongodin et al., 2006). In contrast, in group G5, Methanobrevibacter was found to be more abundant. Methanobrevibacter, a well-known genus within the phylum Euryarchaeota, is a dominant archaeal group in the ruminant gut and plays a crucial role in methane production, affecting the host’s energy metabolism and contributing to greenhouse gas emissions
(Henderson et al., 2015; Janssen and Kirs, 2008). This finding underscores the need for strategies to manipulate gut microbiota to mitigate methane emissions from ruminants. However, the varying abundance of Micrococcaceae;g__Arthrobacter, Methanobacteriaceae; g__Methanobrevibacter, Clostridiaceae;g__Clostridium and Bacteroidaceae; g__Bacteroides across different samples also highlight the intricate, dynamic and personalized nature of the microbiome
(Lloyd-Price et al., 2017). This study provides a detailed snapshot of the taxonomic composition of the yak microbiome, but further studies involving larger sample sizes and various environmental conditions are required to fully understand the factors that govern this diversity.
Kyoto encyclopedia of genes and genomes (KEGG), evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) and Carbohydrate-Active enZymes (CAZy) databases provide invaluable resources for the functional annotation of metagenomic sequences, which contribute significantly to the understanding of the functional potential and metabolic capacity of microbiomes. The functional analysis of unigenes in the yak microbiota using the KEGG database reveals an intriguing picture of the microbial functions in the gastrointestinal system of these animals. A high number of unigenes were linked to cellular processes, environmental information processing, genetic information processing, human diseases, metabolism and organismal systems. These findings reflect the profound influence of gut microbiota on host health and physiology, substantiated by earlier research
(Nicholson et al., 2012; Turnbaugh et al., 2006). The KEGG analysis delineated noticeable differences among the five study groups, particularly group G4, which displayed more metabolic genes and fewer disease-related genes (Fig 9). This result corroborates the widely held view that a balanced gut microbiota with high metabolic capability can contribute to host health and may confer protection against certain diseases (
Honda and Littman, 2016). G4’s increased metabolic functions could be driven by several factors, including diet, host genetics, or environmental factors, more energy for muscle activity associated with male, all of which have been shown to influence the microbial community structure and function
(Lozupone et al., 2012). Our statistical analysis, including PCA, PoCA and ANOSIM, further emphasized the functional differences among the groups (Fig 10). These results align with previous work demonstrating substantial inter-individual and intra-group variability in microbiome functions
(Falony et al., 2016). However, the specific factors leading to these differences in our study are not clear and would require further investigation. The minor differences in our study compared to some other published works may be due to several reasons, such as variations in sampling and data analysis methods, differences in host genetics and environmental factors
(Arumugam et al., 2011). Moreover, the specificity of yak gut microbiota, shaped by unique dietary and environmental pressures, could also lead to functional profiles that differ from those found in other animals or human studies
(Ley et al., 2008).
Further the present study also conducted an analysis of unigenes from the yak microbiota using the eggNOG and CAZy databases that provides an insightful view of the functional attributes of the yak microbiota, focusing on various aspects of metabolism and carbohydrate processing. The eggNOG database revealed a significant presence of genes linked to various metabolic pathways such as energy production and conversion, amino acid transport and metabolism, nucleotide transport and metabolism, among others. This aligns with previous research, which indicates that gut microbiota plays a crucial role in the host’s metabolic processes (
Bäckhed et al., 2005). Notably, Group G4 exhibited higher abundance of amino acid transport and metabolic genes. It’s well established that gut microbes can contribute to the host’s amino acid metabolism
(Neis et al., 2015). However, the precise factors driving the higher abundance of these genes in Group G4 warrant further investigation. Our analysis with the CAZy database showed a significant number of unigenes involved in various carbohydrate processing functions such as glycoside hydrolases, carbohydrate”binding modules and carbohydrate esterases. Glycoside hydrolases, especially, play an important role in the breakdown of dietary fibers in the gut
(Flint et al., 2008). Interestingly, group 4 had fewer glycoside hydrolases genes, while group G3 had the most. This may suggest differences in dietary fiber processing abilities between these groups. Further, these findings could also indicate differences in diet or host-specific factors that affect the composition of the microbiota
(David et al., 2014). Our PCA, PoCA and ANOSIM analysis demonstrated functional and metabolic differences between the groups, resonating with prior studies indicating considerable variability in microbial functions across individuals and groups
(Falony et al., 2016). Comparatively, some differences between this and other studies could be due to diverse factors such as methodology, host genetics, environmental factors and the specific diet and lifestyle of yaks. Future investigations might provide more context to these results and contribute to a better understanding of the yak microbiota.
The study of antibiotic resistance genes (ARGs) in the microbiome of yaks of different ages bears significant implications for public and veterinary health and for the broader understanding of the resistome - the collective reservoir of antibiotic resistance in microbial communities. It permits the tracking of the development and spread of antibiotic resistance over time and changes in ARG composition can serve as a bellwether for future trends in resistance (
O’Toole and Jeffery, 2015;
Van Boeckel et al., 2015). Furthermore, this research can guide the practice of antibiotic stewardship in veterinary medicine, helping to slow the advent of further resistance
(Lhermie et al., 2017). we investigating the distribution of antibiotic resistance genes (ARGs) within the microbiome of yaks at various ages, our findings illustrate the dynamic and complex nature of ARGs within microbial communities. In G1, Bacillota accounted for 79% of ARGs, while Bacteroidota and Actinomycetota accounted for 5% and 3%, respectively. The relative abundance analysis results in the present study suggest that there are clear differences in the antibiotic resistance ontology (ARO) among different yak groups. In terms of the top 20 AROs, the relative abundance of vancomycin resistance genes was found to be the highest in all groups, followed by resistance genes for macrolides, tetracycline, oxazolidinones, bacitracin and aminoglycosides. It is interesting to note that the vanW_gene_in_vanI_cluster was found to be highest in the G1 group and lowest in the G4 group. This difference could be influenced by several factors including environmental factors that could influence ARGs
(Shenhav et al., 2019). In contrast, the G4 group showed the highest relative abundance of vanY_gene_in_vanB_cluster, suggesting differential antibiotic pressures or environmental exposures for this group. Heat map results and ANOSIM analysis based on Bray-Curtis distance further underline these area-related differences. The clear clustering of different area groups in the heat map reflects distinct differences in their ARG profiles. ANOSIM analysis, a non-parametric test used to compare the similarity of microbial communities, further confirms these findings, showing that the between-group differences are significantly higher than the within-group differences. These results underline the complexity and heterogeneity of ARGs in different groups of yaks and highlight the influence of multiple factors, including age, sex and possibly antibiotic exposure in specific areas (
Kim et al., 2011).
The dissemination of antibiotic resistance genes (ARGs) in the environment and among different species is indeed a critical public health concern. Birds, as well as companion animals, can act as vectors for the transmission of ARGs due to their ability to migrate or through close contact with humans and other animals, respectively (
Allen et al., 2010;
Cormier et al., 2019). Environmental sources like soil, rivers and sediments serve as reservoirs for ARGs. The horizontal gene transfer, mainly via mobile genetic elements like plasmids, allows ARGs to move between bacteria in these environments. Therefore, animals living in such environments can acquire these ARGs, even in the absence of direct antibiotic exposure
(Bengtsson-Palme et al., 2018). Yaks, being grazing animals, continuously interact with the soil ecosystem while foraging. They can thereby ingest soil bacteria carrying ARGs, which then become incorporated into their gut microbiota. Moreover, as yaks drink water from rivers and other freshwater sources, they could be ingesting water-borne bacteria containing ARGs. These environmental influences could explain why the intestinal microbiome of yaks harbors ARGs despite not being directly treated with antibiotics. The finding underscores the importance of holistic, One Health approaches to tackling antibiotic resistance that consider not just medical and agricultural antibiotic use, but also environmental reservoirs of resistance.
The gastrointestinal tract of animals, including yaks, harbors a complex ecosystem of microorganisms known as the gut microbiota. These microorganisms play a crucial role in host health, including nutrient metabolism, immune system development and protection against pathogens. Understanding the composition and functional potentials of the gut microbiota in yaks is essential for gaining insights into their digestive physiology, overall well-being and performance. A metagenomic sequencing approach was employed to investigate the gut microbial diversity, drug resistance, composition and functional potentials of the gut microbiota. The results from the metagenomic analysis showed a diverse range of microorganisms present in the yak gut microbiota that varied among different yak age groups. The predominant phyla were Bacillota, Bacteroidetes, Proteobacteria and Actinobacteria. The KEGG metabolic pathway prediction shows that amino acid metabolism and carbohydrate metabolism were abundant. The study also revealed an alarming abundance of antibiotic-resistance genes, suggesting that the gut microbiota of yaks could potentially serve as a reservoir for antibiotic resistance, which has significant implications for public health. By applying this approach to fecal samples from yaks of different age groups, this study aimed to unravel the microbial diversity, functional potential and drug resistance in the yak gut microbiome. The identified targets can be used as prebiotics and/or probiotics to improve the overall well-being and production by manipulating the dairy feed with desired microbes
(El-Hentati et al., 2023; Gangil et al., 2021; Maftei et al., 2022). 
