Bayesian modeling reveals host genetics associated with rumen microbiota jointly influence methane emission in dairy cows

The study addresses how host genetics and rumen microbiota jointly influence methane emissions from dairy cows, a major contributor to greenhouse gases and a loss of dietary energy. Prior work has shown methane production is shaped by ruminal microbiota, diet, physiology, and host genome, with archaea central to methanogenesis. Earlier studies suggested genetic influence on methane yield and microbiota, but robust biomarkers have been elusive due to limited sample sizes and inadequate statistical models. The authors aim to quantify the contributions of host genetics and microbiota to methane emission variation and identify host genetic factors affecting microbiota composition. They hypothesize: (1) a Bayesian model allowing joint effects of host and all bacterial genera can quantify microbiota’s contribution to methane emission; (2) host genetics and specific host genes influence microbiota composition; (3) microbial abundances influence methane emission; and (4) both host genetics and microbiota contribute to methane emission variation.

Previous research has linked ruminal microbial composition, feed intake, diet, physiology, and host genome with methane emission. Roehe et al. identified microbial genes associated with methane yield and suggested the archaea:bacteria ratio may be heritable and linked to methane yield, indicating potential for selection based on microbiota. Ross et al. and Difford et al. modeled microbiota effects using similarity matrices as random effects in linear mixed models. However, those approaches implicitly assume equal contribution of microbial components. A Bayesian mixture framework can jointly model individual microbial effects with different shrinkage levels, better capturing the complex structure and sparsity of microbiota effects.

Ethics: All procedures followed the approved protocol (Approval number 2016–15–0201–00959, Danish authorities). Participants: 750 Holstein cows with methane phenotypes; 691 genotyped on Illumina BovineSNP50 v2. Genotype QC: Excluded SNPs with unknown/non-autosomal location, individual call rate <99%, SNP call rate <98%, Hardy–Weinberg p<1e-8, MAF<0.02, retaining 39,034 autosomal SNPs. Pedigree traced back to 1926. Methane phenotyping: CH4 and CO2 concentrations measured spectroscopically during automated milking; means corrected for systematic effects (diurnal, day-to-day). Total methane emission derived using CH4:CO2 ratio and predicted CO2 production from heat production units. Microbiota sampling and processing: Rumen liquid collected via oral rumen scoop. DNA extracted (Qiagen QIAamp). 16S rRNA sequencing with bacterial V1–V3 and archaeal V4–V6 primers; 96 libraries sequenced (Illumina MiSeq and HiSeq, 250 bp paired-end). Reads processed with quality control; clustered into OTUs using Lotus. Consensus sequences built; phylogenetic tree constructed. Taxonomy assigned via RDP classifier (confidence 0.8) using Greengenes gg_13_8_otus. Retained OTUs classified as k_Bacteria or k_Archaea. Data transformations for OTUs: Kept OTUs present in ≥50% of cows; added 0.001 to counts, natural log-transformed, standardized to mean 0, variance 1. To address compositionality, also applied centered log-ratio (CLR) transformation using CoDaSeq; both transformations analyzed separately. Bayesian four-component mixture model: Jointly modeled methane emission as a function of fixed effects, host polygenic effect, and OTU effects. Response y: methane production (l/day). Fixed covariates X: intercept, herd (6 levels), parity (4), days in milk (1–350), and Wilmink function terms for early lactation nonlinearity. Polygenic effect u ~ N(0, Aσu²), where A is pedigree-based additive relationship matrix. OTU matrix M contains transformed abundances. OTU effects β follow a four-component mixture prior (tiny/small/medium/large effects) with mixture proportions (π1..π4) ~ Dirichlet(125,25,5,1) to induce sparsity. Variance components for mixture classes impose increasing variances (strong-to-weak shrinkage). Residuals assumed iid Normal with variance σe². Posterior inference via MCMC (Bayz): 100,000 iterations, 10,000 burn-in, thinning every 20 iterations (4,500 posterior samples). Computed var(Zu) and var(Mβ) each iteration; summarized by posterior means and 95% highest posterior density (HPD) intervals. Mixed model REML for heritable microbiota: Estimated genomic heritability (variance explained by SNPs) for microbial traits (genus-level abundances and selected OTUs), as well as diversity summaries. Built genomic relationship matrix (G) from 50k SNPs per Yang et al. Model: y = 1μ + Zu + e, with u ~ N(0, Gσu²), e ~ N(0, Iσe²); h² = σu²/(σu²+σe²). Tested significance using REML SE; multiple testing controlled using Benjamini–Yekutieli FDR (adjusted p<0.05). Microbiota diversity and composition summaries: Bray–Curtis PCoA and Chao1 to assess beta and alpha diversity on log-transformed standardized data; PCA on CLR-transformed data with Euclidean distances. SNP–microbiota association (GWAS-style mixed model): y = 1μ + Zu + Xg + e, with u ~ N(0, Gσu²); tested each SNP (allele dosage 0/1/2) for association with microbial traits (PCoA axes, Chao1, genera, selected OTUs, and CLR PCs). p-values from t-tests on SNP effects; Bonferroni correction applied across markers.

• Bayesian joint modeling (log-standardized OTUs): Host genetics explained 22% of methane variation (95% HPD: 3–45%). Microbial OTUs explained 7% (95% HPD: 0.02–17%). When fitted jointly, total explained variance was 31%, with host genetics 24% (95% HPD: 3–48%) and microbiota 7% (95% HPD: 0.06–17%).
• Bayesian joint modeling (CLR-transformed OTUs): Separately, host genetics explained 22% (95% HPD: 4–42%) and microbiota 11% (95% HPD: 0.5–22%). Jointly, host genetics explained 11% (95% HPD: 0.9–23%) and microbiota 24% (95% HPD: 5–46%).
• Contribution of individual taxa: Each genus/OTU explained a small portion of methane variance; the largest single effect was a bacterial genus (Actinobacteria; Coriobacteriia; Coriobacteriales; Coriobacteriaceae) explaining 1.6% (95% HPD: 0.004–2.4%) of total variance.
• Heritability of microbiota (REML): Highest point estimates observed were h² ≈ 33% for a bacterial OTU, 26% for an archaeal OTU, and 41% for a microbial genus. After multiple testing correction (Benjamini–Yekutieli), none remained statistically significant.
• Host SNP associations with microbiota: Mixed model SNP tests identified markers within host genes ABS4 and DNAJC10 associated with microbiota composition (surviving multiple-testing correction as reported for SNP-level tests).

Findings support that both host genome and rumen microbiota contribute meaningfully to between-cow variation in methane emissions. The host genetic component is substantial and comparable to or exceeding the microbiota component depending on data transformation. Despite clear joint contributions, individual microbial taxa each explain only small fractions of methane variance, indicating methane emission is influenced by many microbial features with small effects. The Bayesian mixture approach effectively modeled sparse, high-dimensional OTU effects while accounting for host polygenic effects, revealing that considering compositionality (CLR) can shift variance attribution between host and microbiota. Evidence of host genomic control over microbiota (moderate heritability estimates and SNP associations in ABS4 and DNAJC10) supports a host–microbiome axis affecting methane emission, aligning with the hypothesis that selecting for host genotypes could indirectly shape microbiota toward lower methane phenotypes. Overall, the integrated modeling provides a framework to dissect host–microbiome interactions relevant to climate mitigation and feed efficiency in dairy systems.

This study demonstrates that Bayesian joint modeling can partition methane emission variance into host genetic and microbiota components, with both contributing appreciably. Host genetics explained roughly 11–24% and microbiota 7–24% of variation depending on transformation, while individual taxa had small effects (largest ≈1.6%). Although point estimates suggested some microbial taxa are moderately heritable, none remained significant after multiple-test correction, emphasizing the polygenic and complex nature of microbiome control. Associations of host SNPs within ABS4 and DNAJC10 with microbiota composition indicate specific host loci may influence the rumen ecosystem. These insights suggest breeding strategies targeting host genetics could help establish favorable microbiota associated with reduced methane emissions. Future work should employ larger cohorts, longitudinal sampling, multi-omics (metagenomics/metatranscriptomics), and refined compositional methods to validate loci, improve biomarker robustness, and translate findings into genetic selection tools.

• After correcting for multiple testing across many microbial traits, no heritability estimates remained statistically significant, indicating limited power and potential false positives among point estimates.
• Individual microbial taxa explained very small fractions of methane variation, complicating biomarker discovery and suggesting a highly polygenic, diffuse architecture.
• 16S rRNA-based profiling limits functional resolution and may introduce biases; platform differences (MiSeq vs HiSeq) were present, though addressed in modeling.
• Microbiota data are zero-inflated and compositional; although CLR was analyzed, residual compositional and sparsity issues may affect estimates.
• HPD intervals for variance components were wide in places, reflecting uncertainty.
• Cross-sectional design limits causal inference; diet and management covariates were adjusted, but unmeasured environmental factors may confound results.