Seven clinically healthy racing dromedary camels (Camelus dromedarius), aged 4–6 years, were selected from a private racing camel facility in Dubai, UAE. The animals were maintained under uniform management and feeding conditions. The diet consisted of a balanced concentrate mix and 6 mm pelleted feed (Pellet I) formulated to meet the nutritional requirements of athletic camels (Table 1). Both feed types supplied sufficient energy and fiber for optimal rumen fermentation and performance.

Rumen contents were collected from each camel before morning feeding using a sterile flexible stomach tube (10). Approximately 500 mL of rumen fluid was obtained and filtered through two layers of sterile cheesecloth to remove large particles. The pH was measured immediately using pH indicator strips (range 4.0–7.0) to confirm rumen health status. Aliquots of the clarified rumen liquor were snap-frozen in liquid nitrogen on-site, transported to the laboratory on dry ice, and stored at −80 °C until DNA extraction.

Total microbial genomic DNA was extracted from 200 mg of rumen content using the QIAamp^® PowerFecal^® DNA Kit (Qiagen, USA), following the manufacturer's instructions (13). The quality and integrity of extracted DNA were assessed using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific) and Agilent Fragment Analyzer 5400 with the Genomic DNA Analysis Kit (DNF-488). DNA concentration was then quantified using a Qubit Fluorometer (Thermo Fisher Scientific) and the Qubit dsDNA HS Assay Kit.

High-quality DNA was used for library preparation with the Novogene NGS DNA Library Prep Set (catalog no.PT004) according to the manufacturer's protocol. Libraries were assessed for quality and quantity using the KAPA Library Quantification Kit (Roche) and Qubit dsDNA HS Assay Kit. Library size distribution was evaluated using the Agilent 2100 Bioanalyzer with the Agilent High Sensitivity DNA Kit, and qualitatively with the Agilent Fragment Analyzer system using the Qualitative DNA Kit (DNF-915-K1000). Sequencing was performed on the Illumina NovaSeq X Plus platform, generating 150 bp paired-end reads (PE150) using a 300-cycle reagent kit.

Alpha diversity metrics including Species richness, Shannon, Simpson, Evenness and Chao1 indices were computed using the scikit-bio package in Python to evaluate within-sample archaeal diversity (22). Beta diversity was calculated using both Bray–Curtis dissimilarity and Jaccard distance to assess between-sample community variation (23). Principal Coordinate Analysis (PCoA) was applied to visualize compositional dissimilarities, and hierarchical clustering heatmaps were generated to depict archaeal distribution patterns across samples (24).

Functional diversity metrics (functional richness, Shannon and Simpson) were computed based on KO abundance profiles. The archaeal core species and core functional features were defined as those present in 80 % and 100 % of samples with a relative abundance threshold ≥0.01%, respectively. The prevalence and abundance of core functions were visualized using heatmaps generated in seaborn (25).

To assess archaeal stability across animals, mean relative abundance and coefficient of variation (CV) were calculated for each archaeal species (26). Taxa with high abundance and low CV were identified as stable core members. For exploratory differential abundance, taxa and gene families with high variation were highlighted for potential functional specialization among individuals.

Descriptive statistical analyses were performed to summarize the archaeal taxonomic and functional metrics across all samples, including relative abundances, alpha diversity indices, beta diversity distances, and functional diversity indices expressed as mean ± SD. All computational analyses and visualizations were conducted in Python using the libraries pandas, numpy, matplotlib, seaborn, and scikit-bio. Heatmaps, prevalence curves, UpSet plots, and stability plots to depict abundance distributions, core community composition, taxa overlap, and inter-sample variability (27). Since all samples represented healthy camels without treatment grouping, no inferential statistical tests were applied; instead, results were summarized descriptively to reflect conserved features in the camel archaeal community.

A total of 14 paired-end shotgun metagenomic libraries were generated from seven camel rumen samples. Sequencing quality metrics from the Nephele platform indicated high-quality datasets with low duplication rates and consistent GC content across samples. The number of raw reads per replicate ranged from 46.9 million to 75.2 million (mean = 64.5 million reads per sample). The average GC content was approximately 48–50%, consistent with rumen-associated archaeal and bacterial genomes. Duplication rates were low, ranging between 16.4% and 20.4%, suggesting minimal amplification bias. Following read preprocessing, adapter trimming, and quality filtering, host-derived reads were removed by aligning sequences to the Camelus dromedarius reference genome. After host read removal and processing through the Nephele bioBakery WGS pipeline, archaeal reads were successfully identified in all samples. Across the seven camel rumen samples, the total number of archaeal reads ranged from 0.78 to 2.47 million per sample, with an average of approximately 1.75 ± 0.55 million reads. Despite differences in total sequencing depth, archaeal populations were consistently detected in all samples, reflecting the stable presence of methanogenic lineages within the rumen ecosystem.

At the phylum level, the archaeal community was remarkably consistent across all seven camel rumen samples (Figure 1A, Supplementary Table S1). The phylum Euryarchaeota dominated with a mean relative abundance of 50.10 ± 0.02 %, while the Methanomada group closely followed at 49.67 ± 0.08 %. All other archaeal phyla Stenosarchaea group (0.22 ± 0.06 %), Nitrososphaerota (0.005 ± 0.003 %), Thermoproteota (0.004 ± 0.002 %), Nitrososphaerota incertae sedis (0.001 ± 0.002 %), and Nanoarchaeota (undetectable or ≈ 0.000 %) collectively accounted for less than 1 % of archaeal reads in any sample. At the family level, the archaeal community in the camel rumen was dominated by the Methanobacteriaceae, which accounted for a mean relative abundance of 99.58 ± 0.09 % across the seven samples. The next most abundant families were Methanomassiliicoccaceae (mean = 0.062 ± 0.04 %) and Methanosarcinaceae (mean = 0.061 ± 0.03 %). All other archaeal families detected had individual mean abundances below ~0.03 % and collectively contributed less than 0.3 % of the total archaeal reads per sample (Figure 1B, Supplementary Table S1). At the genus level, archaeal community composition was strongly dominated by Methanobrevibacter, accounting for an average of 91.19 ± 3.13 % of the archaeal reads across all seven camel rumen samples. The second most abundant genus, Methanosphaera, comprised 8.38 ± 4.30 % of reads. All other archaeal genera individually contributed ≤ 0.12 % of reads and collectively represented less than 2 % of the archaeal community in any sample (Figure 2A, Supplementary Table S1). At the species level, archaeal community in the camel rumen was dominated by Methanobrevibacter sp. YE315 and Methanobrevibacter millerae, which together accounted for approximately 50.64% ± 4.23% of total archaeal reads across samples. Methanobrevibacter sp. YE315 averaged 34.81 ± 4.68% while Methanobrevibacter millerae averaged 16.1 ± 3.13% Secondary species such as Methanobrevibacter ruminantium (12.93 ± 2.09%) and Methanosphaera sp. BMS (8.06 ± 1.11%) were present at lower but consistent levels, while all remaining detected species together comprised less than 6% of the archaeal population in any individual sample (Figure 2B, Supplementary Table S1).

Alpha-diversity analyses of the archaeal communities revealed consistent within-sample diversity across the seven camel rumen samples (Figure 3). Observed species richness ranged from 271 to 380 taxa, with a mean of 335.14 ± 39.1 across the seven camel rumen samples. The Shannon index averaged 2.346 ± 0.024, and the Simpson index was 0.887 ± 0.002 across the seven camel rumen samples. The evenness index averaged 0.404 ± 0.006, indicating a relatively uniform distribution of archaeal taxa, while the Chao1 richness estimator averaged 429.1 ± 68.5, reflecting high species richness across the seven camel rumen samples.

Between-sample community dissimilarities were assessed using PCoA based on Jaccard (presence/absence) and Bray–Curtis (abundance) distances (Figure 4). In the Jaccard ordination (Figure 4A), PC1 and PC2 explained 24.9% and 18.9% of total variance, respectively, and sample points were tightly clustered with no host-specific separation. The Bray–Curtis ordination (Figure 4B) yielded greater explanatory power with PC1 at 56.0% and PC2 at 36.5%; still, most samples overlapped substantially, with only minor divergence of Sample 1 and Sample 6.

Hierarchical clustering of the top 30 archaeal species further illustrated the structural consistency of the archaeal community (Figure 5). The heatmap shows two principal clusters: one dominated by members of the genus Methanobrevibacter (e.g., M. millerae, M. sp. YE315, M. ruminantium) with high relative abundances across all samples, and a second cluster comprising a suite of low-abundance taxa (e.g., Methanosphaera sp., Methanobacterium spp.) present at < 1 % relative abundance.

Core archaeal taxa were defined as those present in ≥80% of camel rumen samples with a relative abundance ≥0.01%. Based on this criterion, seven dominant archaeal species constituted the core archaeome across all samples (Figure 6, Supplementary Table S2). These taxa included Methanobrevibacter millerae, Methanobrevibacter sp. YE315, Methanobrevibacter ruminantium, Methanosphaera sp. BMS, Methanobrevibacter olleyae, Methanobrevibacter ruminantium M1, and Methanobrevibacter smithii. The Upset-style overlap analysis illustrates that a large pool of archaeal species (391 species-level clades) meets the prevalence threshold, emphasizing a common shared archaeal community across hosts, with very few taxa unique to specific sample subsets (Figure 7). A stability plot further assessed each species' mean relative abundance versus its coefficient of variation (CV) across all samples (Figure 8, Supplementary Table S3). The stability plot shows that dominant taxa such as Methanobrevibacter sp. YE315 and Methanobrevibacter millerae cluster at high mean abundance and low coefficient of variation, thus reinforcing their core classification.

Shotgun metagenomic functional profiling revealed a diverse repertoire of archaeal functional genes across all camel rumen samples. Functional annotation identified numerous core metabolic and translational processes, predominantly associated with ribosomal protein synthesis, energy metabolism, and methanogenesis (Figure 9, Supplementary Table S4). The most abundant functional profiles across all samples were ribosomal structural proteins, including large subunit ribosomal proteins L24e, L39e, L40e, and L29, and small subunit ribosomal proteins S10, S11, S12, S27e, and S28e. Genes encoding key enzymes in methanogenesis such as F4₂0-non-reducing hydrogenase iron-sulfur subunits (EC 1.12.99–1.8.98.6) and methyl-coenzyme M reductase subunits (EC 2.8.4.1) were also detected.

To identify the conserved metabolic capacities within the archaeal community of camel rumen samples, we analyzed the core archaeal functions, defined as KEGG Orthologs (KOs) present in 100% of samples. A total of 289 core KOs were detected, indicating a stable set of functional genes maintained across the archaeal population (Figure 10, Supplementary Table S5). The dominant functions were primarily related to translation machinery and energy metabolism, reflecting essential cellular processes necessary for archaeal survival in the rumen environment. Notably, ribosomal protein–encoding genes represented the majority of the core functions, including large subunit ribosomal proteins (K02924, K02927, K02896, K02948, K02950, K02978, and K02979) and small subunit ribosomal proteins (K02977, K02981, K02983, and K02984), suggesting a strong transcriptional and translational stability across individuals. Energy-related enzymes such as F420-non-reducing hydrogenase (K14127), methyl-coenzyme M reductase beta subunit (K00319), and formylmethanofuran dehydrogenase subunit A (K00200) were also present in all samples, indicating the persistence of methanogenesis and hydrogenotrophic energy pathways as core archaeal functions in the camel rumen.

Functional pathway reconstruction of the archaeal metagenome revealed key metabolic processes predominantly associated with methanogenic activity and amino acid biosynthesis (Figure 11). The pathway methyl-coenzyme M oxidation to CO₂ (PWY-5209) exhibited the highest abundance across all samples (mean = 260.8 ± 59.3). Several amino acid biosynthetic pathways were also highly represented including L-isoleucine biosynthesis I (ILEUSYN-PWY) (mean = 18.3 ± 10.5), L-valine biosynthesis (VALSYN-PWY) (mean = 22.13 ± 8.8), and L-lysine biosynthesis VI (PWY-5097) (mean = 8.7 ± 6.3). Purine nucleotide biosynthetic routes such as 5-aminoimidazole ribonucleotide biosynthesis II (PWY-6122) (mean = 25.06 ± 10.2), superpathway of 5-aminoimidazole ribonucleotide biosynthesis (PWY-6277) (mean = 24.18 ± 10.6), guanosine ribonucleotides de novo biosynthesis (PWY-7221) (mean = 18.1 ± 8.4), and inosine5′-phosphate biosynthesis III (PWY-7234) (mean = 7.4 ± 5.5) were consistently detected, reflecting the active nucleotide metabolism of rumen archaea. Additionally, factor 420 biosynthesis II (PWY-5198) (mean = 17.3± 9.9) was identified, a pathway critical for electron transfer in methanogenic archaea.

The functional diversity of archaeal communities in the camel rumen, evaluated via richness, Shannon and Simpson indices, revealed a stable gene-repertoire across all samples (mean ± SE: richness = 507 ± 29, Shannon = 5.03 ± 0.05, Simpson = 0.98 ± 0.002) (Figure 12). The functional beta diversity of archaeal communities across the camel rumen samples was examined using Bray–Curtis dissimilarity metrics derived from KEGG ortholog abundance data. Principal coordinate analysis (Figure 13A) revealed a clear clustering pattern among samples, with the first two axes explaining 50.45% and 21.15% of the total variance, respectively. The ordination plot indicated moderate functional differentiation between samples. The hierarchical clustering heatmap of Bray–Curtis dissimilarities (Figure 13B) further supported these findings, illustrating a generally conserved functional landscape among samples, with only subtle variation in pairwise distances (range = 0.07–0.18).