The study incorporated 48 bovine rumen microbial samples encompassing both metagenomic and meta-transcriptomic sequencing datasets (24). The raw sequencing dataset was obtained from the NCBI Sequence Read Archive (SRA) under accession number PRJNA393057. Raw sequencing data underwent preprocessing through the Kneaddata pipeline (v0.7.2) with three critical phases: Quality trimming and adapter removal were initially performed using Trimmomatic (v0.39) (25), followed by host-derived sequence elimination through alignment against the bovine reference genome (GenBank accession: GCF_002263795) via Bowtie 2 (Version 2.5.3) (26). To address elevated ribosomal RNA (rRNA) abundance in meta-transcriptomic profiles, SortMeRNA (v4.3.2) (27) paired with the SMR v4.3 refined database were used to remove rRNA sequences from both data types, thereby mitigating analytical bias in subsequent selected probiotic taxon expression quantification. The resultant high-fidelity cleaned reads served as the foundation for downstream bioinformatic exploration.

The gene abundance and expression were quantified using salmon (v1.3.2; option l A). Methanogenic pathway activity was quantified using a weighted Z-score approach based on metagenomic sequencing data. Genes with transcript per million (TPM) values > 1 in ≥ 10% of samples were retained to ensure robust expression detection. Core methanogenesis-related genes were identified through KEGG module annotation (M00567, M00357, M00356, M00563) using eggNOG database annotations and eggnog-mapper (28).

A dual normalization strategy was employed: gene expression levels were standardized across samples using Z-score transformation: Z g = E g , s − μ g σ g where E g , s denotes the transcript-per-million (TPM) value of gene g in s sample, and μ g (mean) and σ g (standard deviation) were calculated across all samples. This step minimized technical biases arising from heterogeneous expression scales. Subsequently, inverse variance weights ( w g ) were assigned to prioritize genes with stable expression patterns: w g = 1 σ g 2 + ∈ ( ∈ = 10 − 6 to avoid division by zero ) .

Weighted pathway activity for each sample was then computed by integrating Z-scores and weights across all annotated genes in the methanogenesis modules:

where denotes the set of pathway-associated genes. Finally, cross-sample normalization was performed to ensure comparability:

where μ A and σ A are the mean and standard deviation of raw activity scores across all samples. This approach integrates coordinated expression patterns of pathway-associated genes while enhancing statistical robustness through variance-sensitive weighting.

Taxonomic classification of metagenomic and meta-transcriptomic data was performed using Kraken2 (v2.1.2) (29) with parameter of “--paired.” Clean reads were aligned against archaeal reference genomes from the Genome Taxonomy Database (GTDB release 207) (30) which was pre-processed via the Struo2 pipeline (31) with default parameters prior to analysis.

Taxonomic quantification data of archaeal and bacterial species were processed through the QIIME2 pipeline (2024.5) (32). Sequence reads underwent rarefaction to normalize sampling depth, followed by calculation of relative abundance for archaeal and bacterial communities. To assess archaeal diversity, a rarefied species-level count matrix was re-imported into QIIME2 for alpha diversity (Shannon Index).

Random Forest Model for Identifying Key Microbes Associated with Methanogenic Pathway Activity with the following key parameters: mtry = 3, and ntree = 10,000. Between-group differences in diversity metrics were evaluated using the Kruskal-Wallis nonparametric test, while community structure variation was tested via Analysis of Similarities (ANOSIM). Statistical significance was defined at p ≤ 0.05 for all analyses. Data visualization was implemented with the ggplot2 package in R (v4.3.2).

To investigate the relationship between rumen microbial communities and methanogen expression activity in cattle, paired metagenomic and meta-transcriptomic datasets from 48 cattle in Li et al.’s study were analyzed. Taxonomic profiling using Kraken2 (v2.1.2) and the GTDB (release 207) database revealed incongruence between the relative abundance rankings of microbial species in MG and MT datasets (Figures 1a,b). The top five species by metagenomic relative abundance were Methanobrevibacter sp900314635 (3.14%), Prevotella sp900314935 (2.08%), Prevotella sp900316985 (1.36%), Succiniclasticum sp900315925 (1.29%), and Sodaliphilus sp900320055 (1.20%). In contrast, the metatranscriptomic profiles prioritized Treponema D sp004554075 (1.97%), UBA2810 sp002351705 (1.78%), RUG023 sp900315435 (1.24%), Fibrobacter sp001603905 (1.24%), and Treponema D sp902789325 (1.19%). Notably, most highly abundant taxa were assigned GTDB-specific identifiers, indicative of uncultured microbial lineages, suggesting a substantial reservoir of uncultivated microorganisms in the bovine rumen.

Comparative analysis of species-level alpha diversity (Shannon index) demonstrated significantly higher diversity in MG compared to MT datasets (Wilcoxon rank-sum test, p < 0.05; Figure 1c). Previous studies, such as Peng et al., have identified archaeal taxa (e.g., Methanobrevibacter spp.) with high transcriptional activity in domesticated animals. However, the current findings suggest an overall transcriptional suppression state within the rumen microbial community of cattle.

Archaea, recognized as the primary methanogenic microorganisms, are dominated in the bovine gut by genera such as Methanobrevibacter and Methanocaldococcus according to prior studies. Comparative analysis of these genera revealed markedly higher relative abundances in metagenomic profiles compared to their transcriptional activity in meta-transcriptomic (MT) datasets (p < 0.001; Figures 1d,e). Specifically, Methanobrevibacter exhibited a metagenomic abundance of 5.91% (MG) versus 1.11% (MT; Figure 1d), while Methanocaldococcus showed 0.0019 (for MG) vs. 0.00070% (for MT; Figure 1e). This pronounced disparity underscores a systemic transcriptional suppression of methanogenic archaea within the rumen microbial community, aligning with the observed overall reduction in microbial expression activity.

Based on meta-transcriptomic profiling, this study observed that although the overall functional gene expression of dominant methanogens was suppressed, significantly enhanced activity was detected in key methanogenesis-associated metabolic modules. Specifically, the total expression levels of four critical modules in meta-transcriptomes exhibited statistically higher values (p < 0.0001) than their relative abundance in metagenomes: M00356 [Methyl-coenzyme M reductase, core methanogenesis; TPM (MT) = 528.16 vs. TPM (MG) = 174.71; Figure 1f], M00357 [Tetrahydromethanopterin S-methyltransferase, hydrogenotrophic pathway; TPM (MT) = 1647.40 vs. TPM (MG) = 746.32; Figure 1g], M00563 (Acetyl-CoA decarbonylase synthase, acetoclastic pathway; TPM (MT) = 532.78 vs. TPM (MG) = 160.07; Figure 1h), and M00567 (Coenzyme M biosynthesis, methanogen cofactor synthesis; TPM (MT) = 1940.07 vs. TPM (MG) = 808.02; Figure 1i).

Previous studies have often analyzed methanogen abundance as a proxy for methanogenic potential. However, our analysis revealed that methanogen abundance often appeared lower than, or did not consistently correlate with, the functional expression of methanogenesis pathways. Conversely, the expression of genes comprising key methanogenesis-related metabolic modules demonstrated robust activity. Therefore, we established the Methanogenesis Pathway Expression Activity Index based on the gene expression levels of four key KEGG modules associated with methanogenesis (M00567, M00357, M00356, and M00563). This index was used to identify microbial species associated with the activity of these methanogenic pathways.

Random Forest model was employed to identify microbial species explaining variance in MPEAI (Supplementary Figure S1). Among the top 30 species ranked by feature importance, only five corresponded to currently known cultivated species, while the remaining 25 were represented solely by Genome Taxonomy Database (GTDB) identifiers (indicating they are yet-uncultivated). Pearson correlation analysis revealed that 11 of these species exhibited significant negative correlations with MPEAI. These included UBA1240 sp016285185, Treponema D succinifaciens, and nine other yet-uncultivated Treponema species (Supplementary Figure S2).

Due to the prevalence of potentially taxonomically unresolved (as indicated by GTDB identifiers) and currently uncharacterized species among the uncultivated organisms identified, which precludes functional follow-up, we excluded uncultivated species and repeated the Random Forest analysis. Following the exclusion of uncultivated taxa, the Random Forest model identified the three species classified under the genus Evtepia as the top features positively associated with MPEAI (Figure 2a). Species exhibiting significant negative correlations with MPEAI included: Treponema D bryantii D (Figure 2b; R = − 0.449, p = 0.00137), Treponema D bryantii B (Figure 2c; R = − 0.441, p = 0.0013), and Treponema D bryantii A (Figure 2d; R= − 0.486, p < 0.001), Treponema D succinifaciens (Figure 2e; R = − 0.570, p < 0.001), Treponema D porcinum (Figure 2f; R = − 0.439, p = 0.0018), and Treponema D pectinovorum (Figure 2g; R = − 0.363, p = 0.0111). These results indicate a consistent negative association between members of the genus Treponema and MPEAI, suggesting their potential role as probiotics for mitigating methane production.