In this study, a total of 33,925 MAGs from the gut microbiome of Tibetan antelope (Pantholops hodgsonii) were collected, including 26,607 MAGs from the CNGBdb database (Project ID: CNP0001390) and 7,318 MAGs from the Figshare repository (http://doi.org/10.6084/m9.figshare.29376197).

Additionally, 68 metagenomic sequencing samples corresponding to the 7,318 MAGs were retrieved from the National Center for Biotechnology Information (NCBI) under accession number PRJNA1257558, comprising 26 samples from Blastocystis-present individuals and 42 control samples. Blastocystis presence in the original dataset was determined by PCR amplification of the SSU rRNA gene using primers RD5 and BhRDr, followed by bidirectional Sanger sequencing and subtype confirmation through BLAST and PubMLST analyses. All Blastocystis-positive samples were identified as the subtype ST31 (18, 19).

Taxonomic classification of the MAGs was performed using the classify_wf module of GTDB-Tk (v2.3.2) with the GTDB (release 214) (20). Based on the classification results, MAGs assigned to the domain Archaea were selected for further analysis.

To construct a non-redundant archaeal genome catalog, we first assessed the quality of archaeal MAGs using the predict module of CheckM2 (v1.0.1) (21). MAGs with completeness >50% and contamination < 5% were retained (4).

We then employed the sketch function in Mash (v2.2) (22) to generate reduced representations of archaeal genomes for rapid comparison, followed by the dist function to estimate pairwise distances between genomes. When two MAGs showed a distance of zero, only the one with higher completeness and lower contamination was retained.

Finally, we used the run module of GUNC (v1.0.5) (23) to evaluate potential chimerism in archaeal MAGs, and only those with a clade separation score (CSS) ≤ 0.45 were kept (24). This resulted in the construction of a non-redundant, non-chimeric archaeal genome catalog from the gut microbiota of Tibetan antelopes, which was subsequently used for downstream analyses (25).

Protein-coding sequences (CDSs) were predicted and annotated using Prokka (v1.14.6) (26), with the parameter “–kingdom Archaea” to include archaeal-specific non-fragmented proteins from the UniProtKB database, and “—rfam” to enable the identification of non-coding RNAs.

Functional annotation of predicted protein-coding genes was performed by aligning sequences against the Kyoto Encyclopedia of Genes and Genomes (KEGG, Release 20230401) database (27) using DIAMOND (v2.1.8.162) (28) with the parameters: –outfmt 6 –min-score 60 –query-cover 70 –max-target-seqs 5 –masking 1.

The same DIAMOND workflow was used to annotate CAZymes against the Carbohydrate-Active enZYmes (CAZy, Release 20230726) database (29), with parameters: –outfmt 6 –min-score 60 –query-cover 50 –max-target-seqs 5 –masking 1.

The protein catalog was constructed based on CDSs predicted by Prokka from 463 archaeal MAGs. All protein sequences were clustered using MMseqs2 (v7e2840992948ee89dcc336522dc98a74fe0adf00) (30) with the following parameters: –cov-mode 1 -c 0.8 –kmer-per-seq 80 –min-seq-id 0.5. To evaluate the clustering performance, protein sequences were clustered at different identity thresholds, and the number of unique taxonomic families represented within each cluster was calculated and visualized. To reduce potential contamination, singleton proteins that could not be assigned to any cluster were removed from downstream analyses.

A total of 68 raw metagenomic datasets were first quality-filtered using Fastp (v0.23.0) (31) to remove low-quality reads. Host-derived sequences were removed by aligning reads to the host reference genome (GCA_000400835.1_PHO1.0) using Bowtie2 (v2.5.0) (32).

For relative abundance estimation, 20 million high-quality reads were randomly subsampled from each sample and mapped to both the non-redundant archaeal genome catalog and gene catalog using Bowtie2 with default parameters. Read counts were normalized to transcripts per million (TPM). The relative abundance of each functional category was calculated based on the TPM values of genes assigned to that function.

To identify potential viral sequences within the gut archaea of Tibetan antelope, contigs longer than 5,000 bp from 463 non-redundant archaeal genomes were screened (24, 33). First, the ratio of viral to host genes on each contig was assessed using CheckV (v1.0.1) (34). Contigs containing more than 10 host genes or where host genes outnumbered viral genes by more than fivefold were excluded. Proviral fragments were also identified using CheckV.

To enhance the comprehensiveness of viral detection, three complementary methods were employed: (1) viral gene-enriched regions as determined by CheckV; (2) sequences predicted by DeepVirFinder (v1.0.19) (35) with scores >0.90 and P-values < 0.01; and (3) viral sequences identified by VIBRANT (v1.2.1) (36) using default parameters. Contigs meeting any of these criteria were retained as putative viral sequences.

To remove potential bacterial contamination, candidate viral sequences were screened for bacterial single-copy orthologs using BUSCO (Benchmarking Universal Single-Copy Orthologs) (37) combined with hmmsearch. The BUSCO ratio (number of BUSCOs divided by total gene count) was calculated, and sequences with a ratio ≥5% were discarded.

To characterize viral community structure and assign viral operational taxonomic units (vOTUs) to genus- and family-level taxa, pairwise protein sequence comparisons were performed using BLASTN (v2.13.0) with parameters set to an E-value ≤ 1e-5 and a maximum target sequence number of 99,999. Based on these alignments, the proportion of shared genes and average amino acid identity (AAI) between vOTUs were calculated. Following the approach of Nayfach et al. (38), Markov Cluster Algorithm (MCL) was applied for clustering: at the family level, connections with AAI >20% and shared gene content >10% were retained, with an inflation parameter of 1.2; at the genus level, stricter criteria were applied, requiring AAI > 50% and shared gene content > 20%, with an inflation parameter of 2.

Taxonomic annotation of viral genomes was performed by aligning predicted protein sequences against several reference databases, including the Virus-Host DB (May 2024 release), the crAss-like viral protein database, and other publicly available viral protein datasets (39). Alignments were conducted using DIAMOND with a minimum sequence identity of 30%, query and subject coverage thresholds of at least 50%, and a minimum alignment score of 50. For viral genomes encoding fewer than 30 proteins, assignment to a viral family required at least 20% of the proteins to match proteins from the same family. For genomes encoding 30 or more proteins, at least 10 proteins had to match the same family to allow taxonomic classification (33).

Importantly, none of the viral sequences identified in this study could be assigned to any known viral taxa, suggesting that they may represent novel, previously uncharacterized viral lineages.

CDSs within the identified viral contigs were predicted using Prodigal (v2.6.3) (40) with the parameter “-p meta,” optimized for metagenomic data. The resulting protein sequences were then aligned to the KEGG and CAZy databases using the same annotation strategy and parameters described above, in order to identify the potential metabolic and functional profiles encoded by the viral genomes.

In this study, we compiled a total of 60,664 MAGs from the human gut microbiome (sourced from: https://github.com/snayfach/IGGdb) and 26,053 MAGs from the chicken gut microbiome. The chicken-derived MAGs included 12,339 from the National Microbiology Data Center (NMDC, https://nmdc.cn/icrggc), 979 from the ENA database (project IDs: PRJEB55375 and PRJEB55374), and 6,786 from the Figshare repository (https://dx.doi.org/10.6084/m9.figshare.24901884, https://dx.doi.org/10.6084/m9.figshare.24901878, https://doi.org/10.6084/m9.figshare.24681096.v1). Additionally, 5,949 MAGs were obtained from the NCBI database (project ID: PRJNA1099794).

All human- and chicken-derived MAGs were taxonomically annotated using the same classification pipeline applied to the Tibetan antelope MAGs. Quality filtering and dereplication procedures followed the same criteria used to construct the non-redundant archaeal genome catalog of Tibetan antelope. The resulting archaeal MAGs were then used for comparative analysis of gut archaeal species composition across the three host types.

Phylogenetic trees were visualized using the iTOL platform (https://itol.embl.de) (41). Genome structure and target gene annotations were visualized using the CGView platform (https://proksee.ca/) (42). All statistical analyses were conducted in the R environment (v4.2.1). The heatmap of protein distribution was generated using the “pheatmap” package (v1.0.12) (43). To explore differences in protein composition among archaeal taxonomic groups, non-metric multidimensional scaling (NMDS) was performed for visualization. Based on the presence/absence matrix of protein clusters, Bray-Curtis dissimilarity was calculated using the vegdist function from the “vegan” package (v2.6-4) (44), and NMDS ordination was conducted using the metaMDS function. Stress values were reported for each NMDS plot to evaluate the goodness of fit. The statistical significance of compositional differences among groups was assessed using permutational multivariate analysis of variance (PERMANOVA, 999 permutations) via the adonis2 function in vegan. Rarefaction curves were generated using the R package “vegan.” The Observed index, Chao1 index, and Shannon index were calculated based on relative abundance data. beta-diversity was assessed using Bray-Curtis distances and visualized via principal coordinate analysis (PCoA). Differences between groups were tested using PERMANOVA with 999 permutations to assess statistical significance and quantify the proportion of explained variance. Pairwise comparisons between groups were conducted using the Wilcoxon rank-sum test, with a significance threshold set at P < 0.05. All other visualizations were performed using “ggplot2” (v3.3.6) (45).

A total of 33,925 MAGs were recovered from the gut microbiome of Tibetan antelopes. Among these, 856 archaeal MAGs were identified based on taxonomic annotation. To remove redundancy, we applied Mash-based clustering and obtained 854 non-redundant archaeal MAGs (Figure 1A). We then filtered for medium-quality MAGs with >50% completeness and < 5% contamination, resulting in 567 MAGs. After further removal of chimeric genomes using GUNC, we obtained a final set of 463 non-redundant, non-chimeric archaeal genomes from the Tibetan antelope gut microbiota.

Taxonomically, these MAGs were primarily assigned to the phylum Methanobacteriota (65.22%), followed by Halobacteriota (17.93%) and Thermoplasmatota (16.85%; Figure 1B; Supplementary Table 1). At the genus level, the dominant groups included Methanobrevibacter_A (33.91%), Methanobrevibacter (23.33%), Methanocorpusculum (17.71%), SIG5 (9.72%), Methanosphaera (7.99%), and UBA71 (7.13%). The MAGs exhibited a GC content ranging from 26 to 62% (mean: 39%), and genome sizes from 0.53 Mb to 2.65 Mb (mean: 1.60 Mb). The completeness ranged from 50.55 to 99.97% (mean: 82.30%), while contamination levels varied from 0 to 4.96% (mean: 1.00%).

A total of 691,494 protein sequences were predicted from 463 archaeal MAGs recovered from the gut microbiota of Tibetan antelopes. These proteins were clustered based on >50% amino acid identity and >80% alignment coverage. Singleton clusters containing only one sequence were removed, resulting in 16,189 representative protein clusters. These clusters collectively formed a non-redundant protein catalog for archaeal genomes. Among them, 4,027 clusters were found to be shared across more than 50 MAGs, suggesting the presence of widely distributed proteins among the archaeal genomes. These shared proteins were primarily derived from Methanobacteriaceae, followed by Methanomethylophilaceae and Methanocorpusculaceae, with only a few from Methanosarcinaceae (Figure 2A). To assess differences in protein composition across archaeal phylogenetic groups, we performed non-metric multidimensional scaling (NMDS) analysis. The results showed that the four dominant archaeal families-Methanobacteriaceae, Methanomethylophilaceae, Methanocorpusculaceae, and Methanosarcinaceae were clearly separated in the NMDS space (Figure 2B), indicating notable differences in their protein composition. Given the predominance of Methanobacteriaceae in the archaeal community, we further examined genus-level variation within this family. NMDS analysis revealed a clear separation among the genera Methanobrevibacter_A, Methanobrevibacter, and Methanosphaera (Figure 2C). Notably, although Methanobrevibacter_A and Methanobrevibacter differ only by a suffix in name, they are classified as distinct genera according to the GTDB-Tk taxonomy, which is consistent with the observed differences in their protein composition.

Blastocystis is a protist widely present in the intestines of humans and animals, and its presence may lead to alterations in the structure and function of the host gut microbiota. To determine whether Blastocystis presence affects the diversity and composition of gut archaea in Tibetan antelopes, we compared the archaeal communities between the Blastocystis-present (BP) group and healthy controls (HC) group.

Rarefaction curve analysis indicated that the archaeal richness in both groups approached saturation with increasing sequencing depth (Figure 3A), suggesting adequate coverage of archaeal diversity. Comparisons of alpha diversity indices, including the Observed (Figure 3B), Chao1 (Figure 3C), and Shannon (Figure 3D) indices, showed no significant differences between the two groups (P > 0.05). However, principal coordinates analysis (PCoA) revealed a significant separation in community structure between the BP and HC groups (Figure 3E, R² = 0.0387, P = 0.023).

At the phylum level, the archaeal communities in both groups were dominated by Methanobacteriota (HC: 96.13%, BP: 87.56%), followed by Halobacteriota (HC: 3.80%, BP: 12.26%) and Thermoplasmatota (HC: 0.07%, BP: 0.17%; Supplementary Figure 1A; Supplementary Table 2). Notably, the relative abundance of Methanobacteriota was significantly reduced in the BP group (P < 0.01), while the proportions of Halobacteriota and Thermoplasmatota were significantly increased (P < 0.01; Supplementary Figures 1B, C).

At the genus level, the dominant archaeal taxa in both groups included Methanobrevibacter_A (HC: 51.97%, BP: 35.80%), Methanobrevibacter (HC: 28.12%, BP: 33.25%), and Methanosphaera (HC: 16.04%, BP: 18.51%; Figure 3F). Notably, Methanobrevibacter_A was significantly decreased in the BP group (P < 0.05), whereas Methanocorpusculum (P < 0.05), SIG5 (P < 0.01), and UBA71 (P < 0.01) were significantly enriched following Blastocystis presence (Supplementary Figure 2).

To investigate the functional characteristics of gut archaea in Tibetan antelopes and assess the potential impact of Blastocystis presence, we performed functional annotation and comparative analysis based on the KEGG and CAZy databases.

For KEGG annotation, a total of 445,954 genes were assigned to 2,280 KEGG Orthologs (KOs). Among these, 141,494 genes were associated with Protein families: Genetic Information Processing (B09182), 69,264 with Protein families: Signaling and Cellular Processes (B09183), and 64,451 with Carbohydrate Metabolism (B09101; Figure 4A). In terms of relative abundance, archaeal functions were predominantly related to genetic information processing, signaling, and energy metabolism (Figure 4B). PCoA analysis revealed a significant shift in KEGG-based functional profiles between BP and HC groups (Figure 4C; R² = 0.0445, P = 0.019). Specifically, Among the functional categories with the highest relative abundances, Energy Metabolism (P < 0.01), Nucleotide Metabolism (P < 0.05), and Unclassified: Metabolism (P < 0.01) were significantly reduced in the BP group, while the relative abundance of Translation was significantly increased (P < 0.05; Supplementary Figure 3).

For CAZy annotation, 63,284 genes were assigned to 199 CAZyme families. Of these, 37,676 genes were related to GlycosylTransferases (GT), 17,360 to Glycoside Hydrolases (GH), 4,996 to Carbohydrate-Binding Modules (CBM), 2,217 to Carbohydrate Esterases (CE), 1,026 to Auxiliary Activities (AA), and 9 to Polysaccharide Lyases (PL; Figure 4D). In terms of relative abundance, GTs accounted for more than 50% of the archaeal carbohydrate-active enzymes (CAZymes), followed by GHs and CBMs (Figure 4E). PCoA analysis further indicated a significant alteration in the carbohydrate metabolic potential of archaea between the BP and HC groups (Figure 4F). Specifically, the relative abundances of CBM (P < 0.01), CE (P < 0.05), and PL (P < 0.05) were significantly decreased in the BP group, while the relative abundance of GT was significantly increased (P < 0.05; Supplementary Figure 4).

To investigate the presence of archaeal viruses in the gut microbiota of Tibetan antelopes, we performed virus identification based on the constructed non-redundant archaeal genome catalog. As a result, seven viral sequences were identified, all of which were unclassified viruses. Functional annotation revealed that two viral sequences were assigned to K20276 (bapA), which encodes a large repetitive protein involved in the quorum sensing pathway (ko02024). In addition, five viral sequences were annotated as members of the glycosyltransferase family 4 (GT4), which are commonly associated with the synthesis and modification of polysaccharides.

To explore the compositional characteristics of gut archaea in Tibetan antelopes, we compared the taxonomic distribution of their non-redundant archaeal genomes (MAGs) with those from humans and chickens. At the phylum level, archaeal MAGs from Tibetan antelopes were assigned to three phyla: Methanobacteriota (65.23%), Halobacteriota (17.93%), and Thermoplasmatota (16.85%; Supplementary Figure 5A). In contrast, the majority of archaeal MAGs from humans belonged to Methanobacteriota (81.87%), with much lower proportions of Halobacteriota (1.55%) and Thermoplasmatota (16.58%; Supplementary Figure 5B). In chickens, Halobacteriota was the dominant phylum (52.02%), followed by Thermoplasmatota (22.54%) and Methanobacteriota (25.43%; Supplementary Figure 5C).

At the genus level, the predominant archaeal genera in Tibetan antelopes were Methanobrevibacter_A (33.91%), Methanobrevibacter (23.33%), and Methanocorpusculum (17.71%; Figure 5A). In humans, Methanobrevibacter_A dominated (80.31%), followed by UBA71 (12.44%) and Methanomassiliicoccus_A (2.07%; Figure 5B). In chickens, the dominant genera were Methanocorpusculum (52.02%), UBA71 (20.23%), and Methanobrevibacter_A (19.65%; Figure 5C).

Notably, 71.50% of archaeal MAGs from Tibetan antelopes could not be classified to known species, whereas the proportion of unclassified MAGs at the species level was only 1.55% in humans and 1.73% in chickens (Supplementary Table 4), indicating a higher degree of novelty and unexplored diversity in the gut archaeal community of Tibetan antelopes.