# The fecal microbiota of Holstein cows is heritable and genetically correlated to dairy performances

L. Brulin, S. Ducrocq, J. Estellé, G. Even, S. Martel, S. Merlin, C. Audebert, P. Croiseau, M.-P. Sanchez

► To cite this version:

L. Brulin, S. Ducrocq, J. Estellé, G. Even, S. Martel, et al.. The fecal microbiota of Holstein cows is heritable and genetically correlated to dairy performances. Journal of Dairy Science, 2024, 107 (12), pp.11254-11268. 10.3168/jds.2024-25003. hal-05010432

HAL Id: hal-05010432

https://hal.inrae.fr/hal-05010432v1

Submitted on 28 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The fecal microbiota of Holstein cows is heritable and genetically correlated to dairy performances

L. Brulin,<fnref n="1" />,<fnref n="2" />* S. Ducrocq,<fnref n="1" />,<fnref n="3" /> J. Estellé,<fnref n="2" /> G. Even,<fnref n="1" />,<fnref n="3" /> S. Martel,<fnref n="1" />,<fnref n="3" /> S. Merlin,<fnref n="1" />,<fnref n="3" /> C. Audebert,<fnref n="1" />,<fnref n="3" />

P. Croiseau,<fnref n="2" /> and M.-P. Sanchez<fnref n="2" />

## ABSTRACT

The fecal microbiota of ruminants constitutes a diversified community that has been phenotypically associated with a variety of host phenotypes, such as production and health. To gain a better understanding of the complex and interconnected factors that drive the fecal bacterial community, we have aimed to estimate the genetic parameters of the diversity and composition of the fecal microbiota, including heritabilities, genetic correlations among taxa, and genetic correlations between fecal microbiota features and host phenotypes. To achieve this, we analyzed a large population of 1,875 Holstein cows originating from 144 French commercial herds and routinely recorded for production, SCS, and fertility traits. Fecal samples were collected from the animals and subjected to 16S rRNA gene sequencing, with reads classified into amplicon sequence variants (ASV). The estimated α- and β-diversity indices (i.e., Observed Richness, Shannon index, Bray-Curtis and Jaccard dissimilarity matrices) and the abundances of ASV, genera, families, and phyla, normalized by centered-log ratio, were considered as phenotypes. Genetic parameters were calculated using either univariate or bivariate animal models. Heritabilities estimates, ranging from 0.08 to 0.31 for taxa abundances and β-diversity indices, highlight the influence of the host genetics on the composition of the fecal microbiota. Furthermore, genetic correlations estimated within the microbial community and between microbiota features and host traits reveal the complex networks linking all components of the fecal microbiota together and to their host, thus strengthening the holobiont concept. By estimating the heritabilities of microbiota-associated phenotypes, our study quantifies the impact of the host genetics on the fecal microbiota composition. In addition, genetic correlations between taxonomic groups and

between taxa abundances and host performance suggest potential applications for selective breeding to improve host traits or promote a healthier microbiota.

**Key words:** heritability, genetic correlations, dairy cows

## INTRODUCTION

In the past decades, the microbiotas of the gastrointestinal tract have been largely studied in livestock species, revealing symbiotic relationships with the host and direct effects on a variety of phenotypes (O'Hara et al., 2020; Reyer et al., 2020; Bozzi et al., 2021; Chen et al., 2023). In ruminants, microbial communities play a crucial role in feed digestion. This process is essential because the host's own metabolism alone cannot fully accomplish it, necessitating the action of microorganisms. In cattle, the diversity and abundance of rumen microbiota taxa have shown relevant associations with milk production (Jami et al., 2014; Xue et al., 2018) and feed efficiency (Monteiro et al., 2022; Na and Guan, 2022). Moreover, these microbial communities have been linked to host health in cattle (O'Hara et al., 2020; Xu et al., 2021). Additionally, ruminal microbiota has shown promise in addressing environmental issues, with specific taxa associated with significant variations in methane production (Hill et al., 2016; Difford et al., 2018; Martínez-Álvaro et al., 2022a). Although the rumen microbiota has received substantial attention, there is a growing interest in microbial communities inhabiting the other gastrointestinal compartments. These compartments have unique physiological properties, resulting in organ-specific microbiota and revealing specific interplays with the host. Recent studies have demonstrated phenotypic associations between the diversity and composition of the fecal microbiota and milk production levels in dairy cattle (Monteiro et al., 2022), supporting the informative nature of this microbiota, the sampling of which remains less intrusive.

Additionally, several studies have highlighted host genetic influence on gut microbiota in monogastric animals, as well as in the rumen microbiota in cattle. These

<footnote n="1">¹GD Biotech - Gènes Diffusion, 59000 Lille, France</footnote>

<footnote n="2">²Université Paris-Saclay, INRAE, AgroParisTech, GABI, 78350 Jouy-en-Josas, France</footnote>

<footnote n="3">³PEGASE-Biosciences, Institut Pasteur de Lille, 59019 Lille, France</footnote>

<footnote>Received April 3, 2024.
Accepted August 8, 2024.
*Corresponding author: l.brulin@genesdiffusion.com</footnote>
investigations have unveiled a complex landscape where nonheritable and heritable taxa coexist. Low to moderate heritability estimates have been reported for diversity indices and some taxa abundances in pigs (Bergamaschi et al., 2020; Déru et al., 2022), poultry (Zhou et al., 2023) and cattle rumen (Difford et al., 2018; Li et al., 2019). These findings suggest the bovine fecal microbiota could also be modulated by the host's genetic makeup. Furthermore, its diversity and composition could be genetically correlated with various host traits of interest.

Therefore, in this context, the present study aims to investigate the relationships between the genetics of dairy cows and their fecal microbiota. Leveraging a large cohort of Holstein cows, our objectives included estimating the heritabilities and genetic correlations of diversity indices and taxa abundances, as well as the genetic correlations among fecal microbiota diversity indices, taxa abundances, and phenotypes of interest for the dairy industry.

## MATERIALS AND METHODS

This work was conducted on farm animals reared for commercial purposes in compliance with the French regulation (Code Rural et de la Pêche Maritime) and the European Council Directive 98/58/EC. In accordance with the legislation, no approval from the Institutional Animal Care and Use Committee or ethics committee was necessary, as all the performed sampling operations were part of routine animal manipulations by duly authorized technicians.

### Population

Fecal sampling was conducted on a cohort of 1,875 lactating Holstein cows between May 2020 and November 2022. These cows were distributed across 144 commercial farms in France, with the number of cows sampled per herd ranging from 5 to 33 cows (mean = 13). All sampled animals were lactating cows, with lactation stages ranging from 5 to 715 DIM (mean = 192), and 31% of the cows were in their first lactation. Because the animals were raised in commercial herds, detailed descriptions of specific management practices or rations were limited. However, within each herd, contemporary individuals shared similar farming conditions and diets. Therefore, both the ration and management practices were considered within the herd-sampling effect in the subsequent analyses, as described below in Model [1].

### Sample Collection, DNA Extraction and Sequencing

Approximately 50 g of fecal samples were collected from cows by technicians at Gènes Diffusion (Douai,

France) by applying a digital rectal stimulation with sterile-gloved hands. Samples were conserved using 30 mL sterile tubes containing 15 mL of an ethanol-based preserving liquid at room temperature during transportation. The samples were then stored in a cold chamber (+4°C, 30-40 d) at the Gènes Diffusion research laboratory (Institut Pasteur de Lille, Lille, France) before being processed.

DNA extraction and 16S rRNA gene sequencing were performed on the samples that were previously grouped into 18 laboratory batches. The DNA extraction was performed using the Nucleospin 96 Soil kit (Macherey Nagel, Düren, Germany), following the manufacturer's recommendations. Then, the sequencing library was constructed using a dual-indexed paired-end sequencing strategy that targeted the V3-V4 variable regions of the 16S rRNA gene. To achieve this, 2 PCR steps were applied to the samples. The first PCR used the 5'-Tag sequences 5'-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG-3' and 5'-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG-3' as the forward and reverse primers, respectively, along with the 16S rRNA gene specific sequences that included forward 5'-CCTACGGGNGGCWGCAG-3' and reverse 5'-GACTACHVGGGTATCTAATCC-3' primers. In addition to the Tag sequences, the PCR2 primers contained a locus to index the samples (barcode sequence) and a locus sequence adapter suitable for the Illumina sequencing technology. The sequencing library was paired-end sequenced on the MiSeq platform (Illumina) using the MiSeq Reagent Kit v3, which performed 600 sequencing cycles.

### Bioinformatics

The raw fastq files underwent quality control using FastQC (Andrews, 2010) at the end of each sequencing run. Sequence analyses were performed with the dada2 v1.24-0 R package (Callahan et al., 2016), following the author's recommendations. The reads from each sequencing batch were trimmed based on the Phred score Q20 as a quality threshold. The DADA2 method (Callahan et al., 2016) was employed to group sequences into amplicon sequence variants (ASV). Subsequently, count tables from each sequencing run were merged to generate a consolidated table. The taxonomic assignment of each ASV was then performed in DADA2 with the SILVA v138.1 reference annotation (Quast et al., 2013). The global composition of the samples was then analyzed using the vegan v.2.6-4 (Oksanen et al., 2022) and phyloseq v.1.44.0 (McMurdie and Holmes, 2013) R packages. Additionally, we investigated the composition of the fecal core microbiota of Holstein cows that included all ASV, genera, families, and phyla observed in at least 90% of the animals.
<table><caption>**Table 1.** Summary statistics of the YD of the phenotypes of the Holstein cows that were sampled in the study</caption><thead><tr><th>Phenotype</th><th>N<sup>1</sup></th><th>Mean</th><th>SD</th><th>Minimum</th><th>Maximum</th></tr></thead><tbody><tr><td>MY</td><td>1,704</td><td>14.75</td><td>71.71</td><td>-249.25</td><td>277.07</td></tr><tr><td>FY</td><td>1,704</td><td>71.90</td><td>307.87</td><td>-968.02</td><td>1,266.98</td></tr><tr><td>PY</td><td>1,704</td><td>39.95</td><td>199.48</td><td>-730.58</td><td>779.82</td></tr><tr><td>FC</td><td>1,705</td><td>1.40</td><td>37.97</td><td>-109.32</td><td>151.27</td></tr><tr><td>PC</td><td>1,704</td><td>-0.21</td><td>17.39</td><td>-64.39</td><td>60.49</td></tr><tr><td>SCS</td><td>1,704</td><td>-1.74</td><td>7.31</td><td>-22.83</td><td>26.93</td></tr><tr><td>CCR</td><td>1,746</td><td>1.15</td><td>2.73</td><td>-7.51</td><td>6.91</td></tr><tr><td>CAI<sub>f</sub></td><td>1,750</td><td>-15.88</td><td>205.01</td><td>-940.83</td><td>930.04</td></tr><tr><td>IFF</td><td>1,732</td><td>-63.43</td><td>419.14</td><td>-871.80</td><td>1,807.56</td></tr></tbody></table>

## Definition of Microbiota-Based Phenotypes

To explore the fecal microbiota, 2 categories of microbiota-based phenotypes were defined. First, abundances of ASV, genera, families, and phyla with a prevalence exceeding 60% were submitted to centered-log ratios (CLR) transformation. The CLR transformation was implemented with the transform() function of the microbiome (v.1.24.0) R package (Lahti and Shetty, 2017), which incorporates a pseudo count equal to half the minimum relative abundance. Outliers, defined as values beyond the range of mean ± 4 SD, were removed. Second, diversity indices were estimated from raw counts in a data set that was rarefied to 5,000 reads per animal. At the α-diversity level, the observed richness and the Shannon index were estimated using the phyloseq (v.1.44.0) and vegan (v.2.6-4) R packages, respectively. The β-diversity was estimated using both the Jaccard and Bray-Curtis dissimilarity matrices, obtained with the vegan R package. For each β-diversity index, a nonmetric multidimensional scaling (NMDS) was constructed using the associated dissimilarity matrix. Subsequent analyses focused on the resulting values of the first 2 axes: NMDS1_Bray-Curtis and NMDS2_Bray-Curtis for the Bray-Curtis dissimilarity matrix, and NMDS1_Jaccard and NMDS2_Jaccard for the Jaccard dissimilarity matrix. Outliers, defined as values beyond the range of mean ± 4 SD, were removed in the subsequent analyses.

## Definition of Cow Performance

For each sampled animal, additional phenotypic information was obtained for 9 diverse traits related to various aspects of cows' performance. Five traits were related to milk production and composition: milk yield (MY), fat yield (FY), protein yield (PY), fat content (FC), and protein content (PC). Fertility was assessed through 3 female reproductive traits: lactating cow's conception rate (CCR) which represents the success or failure of each AI, and intervals between calving and first AI (CAI<sub>f</sub>), and between first AI and fertilizing AI (IFF), measuring

the cow's ability to initiate the postpartum cycle and to be fertilized, respectively. Udder health was evaluated using SCS, computed as the mean of log-transformed SCC.

Monthly records during lactation were considered for MY, FY, PY, FC, PC, and SCS, while fertility traits were measured at each lactation. To mitigate the impact of environmental effects, all the analyzed phenotypes were represented as yield deviations (YD). These YD, adjusted for nongenetic effects and averaged per cow for repeated records, were calculated by the French national genetic evaluation systems for the whole Holstein population using the models described at https://interbull.org/ib/geforms. Descriptive statistics of the YD corresponding to the 10 traits are presented in Table 1.

## Estimations of (Co)variance Components and Genetic Parameters

Variance components and genetic parameters were estimated using univariate and bivariate animal models and algorithms implemented in the Wombat software (Meyer, 2007).

For univariate analyses, the following linear animal model was applied:

$$y = X\beta + Z\alpha + e, \quad [1]$$

where y represents the vector of phenotypes for a given trait, including normalized abundance of ASV, genus, family, and phylum, diversity indices, YD for milk production, milk composition, fertility, and udder health; β is the vector of fixed effects; X and Z are the incidence matrices; α ~ N(0, Aσ<sup>2</sup><sub>α</sub>) is the animal genetic effect with A representing the relationship matrix among all individuals calculated from the pedigree (spanning 4 generations and 8,668 individuals), e ~ N(0, Iσ<sup>2</sup><sub>e</sub>) is the vector of random residual effects; σ<sup>2</sup><sub>α</sub> and σ<sup>2</sup><sub>e</sub> are the additive genetic and residual variances, respectively.

For microbiota-based phenotypes (relative abundances and diversity indexes), fixed effects included in

<footnote n="1">¹N = number of cows.</footnote>
the model were determined through prior permutational multivariate ANOVA (PERMANOVA) analyses using the vegan R package. These fixed effects encompassed sampling within the herd (167 levels, with a minimum of 5 cows per herd and sampling session), parity (4 levels: first, second, third, and ≥ fourth parity), lactation stage (24 classes: 15-d stages from 5 DIM to 320 DIM, 45-d stages from 320 DIM to 365 DIM, 40-d stage from 365 DIM to 405 DIM, and ≥405 DIM), and the number of reads per sample. For nonmicrobiota-based phenotypes, already corrected for environmental effects estimated in the entire French Holstein population, the model only included the qualitative fixed effect of birth year to account for the potential genetic trends in these traits.

The heritability of each trait ($h^2$) was calculated as $\left( h^2 = \frac{\hat{\alpha}_{\alpha}^2}{\hat{\alpha}_{\alpha}^2 + \hat{\alpha}_e^2} \right)$ from the estimates of additive genetic ($\hat{\alpha}_{\alpha}^2$) and residual ($\hat{\alpha}_e^2$) variances.

Covariances were estimated using bivariate models, with genetic and residual variances represented as $2 \times 2$ matrices. The additive genetic correlations ($r_{\alpha}$) between microbiota-based phenotypes and between microbiota-based phenotypes and nonmicrobiota-based phenotypes, were calculated as $\left( r_{\alpha} = \frac{\hat{\sigma}_{\alpha1,\alpha2}}{\hat{\sigma}_{\alpha1} \times \hat{\sigma}_{\alpha2}} \right)$, using the estimates of additive genetic covariance between trait 1 and trait 2 ($\hat{\sigma}_{\alpha1,\alpha2}^2$) and the additive genetic standard deviation estimates for trait 1 ($\hat{\sigma}_{\alpha1}$) and trait 2 ($\hat{\sigma}_{\alpha2}$).

We estimated the P-value ($P$) by comparing the test statistics to the Student's t-distribution with the pt() R function (R Core Team, 2022). A heritability or genetic correlation estimate was considered significant if $P \le 0.05$. Due to the large number of independent tests conducted, we further adjusted the P-values using the false discovery rate (FDR) correction method. We computed FDR-corrected q-values for all estimated genetic parameters, that is, heritabilities and genetic correlations, applying an FDR threshold of 10% (FDR 10%).

# RESULTS

## Taxonomic Characterization of the Fecal Microbiota

The analyses of 16S rRNA sequences revealed a rich microbial composition within the fecal microbiota of our population of Holstein cows (n = 1,875), identifying different 37,138 ASV with an average of 19,891 paired-end reads per animal (Supplemental Figure S1; see Notes). These ASV were categorized into 474 genera, 208 families, and 33 phyla. The most abundant ASV (2.0% of total abundance) was annotated as *Paeniclostridium* at the genus level. *Firmicutes* emerged as the predominant

phylum, accounting for 54.9% of the reads (Figure 1A), followed by *Bacteroidetes* as the second most represented phylum, constituting 35.2% of the reads. Among other taxonomic ranks, an unclassified genus from the *UCG-010* family, *Rikenellaceae RC9 gut group*, *Bacteroides*, and *UCG-005* were the most abundant genera, with average relative abundances (RA) of 8.2%, 8.1%, 6.0%, and 6.0%, respectively. The 3 most abundant families were *Rikenellaceae*, *Lachnospiraceae*, and *Oscillospiraceae*, with average RA of 13.6%, 9.8%, and 8.8%, respectively.

We further investigated the composition of the fecal core microbiota of Holstein cows (prevalence ≥60%) and identified 32 ASV, 52 genera, 36 families, and 8 phyla, representing 0.09% of ASV, 11.1% of genera, 17.8% of families, and 24.2% of phyla, respectively (Figure 1B). Interestingly, among these taxa, 5 genera (*Phascolarctobacterium*, *Rikenellaceae RC9 gut group*, an unknown genus from *Lachnospiraceae* family, *Prevotellaceae UCG.003*, and an unknown genus from the order *Oscillospirales*), 7 families (*Acidaminococcaceae*, *Oscillospiraceae*, *Rikenellaceae*, *Lachnospiraceae*, *Prevotellaceae*, *Ruminococcaceae*, and an unknown family from the order *Oscillospirales*), and 3 phyla (*Proteobacteria*, *Bacteroidota*, *Firmicutes*) were consistently observed in all samples. Conversely, no ASV was found to be shared among all samples. The most prevalent ASV, with a 16S rRNA gene closely related to *Romboutsia ilealis* and *Romboutsia timonensis*, exhibited the highest prevalence, reaching 99.5%.

## Heritabilities of Microbiota-Based Phenotypes

Our analyses of microbiota diversity measurements (Table 2) revealed negligible heritability estimates for the $\alpha$-diversity index (i.e., Observed Richness and the Shannon index). However, a significant heritability estimate of 0.14 (SE = 0.06, $P$ = 0.010) was obtained for the second NMDS axis of the Bray-Curtis dissimilarity matrix.

<table><caption>**Table 2.** Heritability estimates for $\alpha$- and $\beta$-diversity indices and their respective SE<fnref n="1" />)</caption><thead><tr><th>Diversity index</th><th>h<sup>2</sup></th><th>SE</th></tr></thead><tbody><tr><td><strong>α-Diversity</strong></td><td></td><td></td></tr><tr><td>Observed richness</td><td>0.028</td><td>NA<sup>2</sup></td></tr><tr><td>Shannon index</td><td>0.011</td><td>0.042</td></tr><tr><td><strong>β-Diversity</strong></td><td></td><td></td></tr><tr><td>NMDS1_Jaccard</td><td>0.064</td><td>0.049</td></tr><tr><td>NMDS2_Jaccard</td><td>0.049</td><td>0.052</td></tr><tr><td>NMDS1_Bray_Curtis</td><td>0.047</td><td>0.045</td></tr><tr><td>NMDS2_Bray_Curtis</td><td><strong>0.135</strong></td><td><strong>0.061</strong></td></tr></tbody></table>

<footnote n="1">¹Bold values indicate significant heritability (FDR 10%).</footnote>

<footnote n="2">²NA = not available due to convergence issues.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Taxonomic characterization of the fecal samples of Holstein cows. (A) Relative abundance of the 10 most abundant phyla of the fecal microbiota. The “Others” group gathers all the other phyla. (B) Proportion of taxa shared by the samples, with the red curve representing the ASV, the green curve for the genus, the blue curve for family, and the purple for the phylum taxonomic ranks. The vertical line highlights the 90% prevalence threshold, used to define the core microbiota.</figcaption></figure>

For estimates on taxa abundances, after filtering for prevalence superior to 60%, we considered 162 ASV, 90 genera, 51 families, and 11 phyla, which were observed in more than 60% of the samples. This represented 0.4% of ASV, 26.9% of genera, 35.9% of families, and 33.3% of phyla. Among them, 15 ASV, 14 genera, 11 families, and 3 phyla presented significant heritability ($P \le 0.05$; Figure 2). However, after applying an FDR 10% correction, only one ASV from the *Negativibacillus* genus remained significant with a heritability estimate of 0.31 (SE = ±0.09). At the genus level, 4 genera were found heritable after correction: *Parasutterella* ($h^2$ = 0.21, SE = 0.07), *Akkermansia* ($h^2$ = 0.21, SE = 0.07), *Paenilostridium* ($h^2$ = 0.21, SE = 0.07), and an unclassified genus from the *WCHB1-41* order ($h^2$ = 0.19, SE = 0.07). At the family level, 5 families presented significant q-values: *Sutterellaceae* exhibited the highest heritability ($h^2$ = 0.25, SE = 0.08), followed by *Akkermansiaceae* ($h^2$ = 0.22, SE = 0.08), an unclassified family from the *WCHB1-41* order ($h^2$ = 0.19, SE = 0.07), *Peptotreptococcaceae* ($h^2$ = 0.18, SE = 0.07), and *Paludibacteraceae* ($h^2$ = 0.17, SE = 0.07).

## Heritabilities of Cow Performance

The heritability estimates of cow performance were highest for production traits. Among these, PC ($h^2$ = 0.73, SE = 0.09, *P* < 0.001) and FC ($h^2$ = 0.72, SE = 0.09, *P* < 0.001) presented the highest values, followed by FY ($h^2$ = 0.60, SE = 0.09, *P* < 0.001), PY ($h^2$ = 0.45, SE = 0.08, *P* < 0.001), and MY ($h^2$ = 0.38, SE = 0.08, *P* < 0.001). Somatic cell score showed a heritability value of 0.43 (SE = 0.08, *P* < 0.001). Finally, the heritability estimates of fertility phenotypes, that is, CAI<sub>f</sub> ($h^2$ = 0.12, SE = 0.06, *P* = 0.023) and CCR ($h^2$ = 0.25, SE = 0.08, *P* = 0.047), were lower. Specifically, IFF ($h^2$ = 0.04, SE = 0.05, *P* = 0.21) was not significantly heritable.

## Genetic Correlations Among Microbiota-Based Phenotypes

We investigated the genetic correlations between heritable taxa within the fecal microbiota (Figure 3). Although, some taxa were uncorrelated within all taxonomic ranks (i.e., ASV, genus, family, and phylum),
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Heritability estimates of the CLR abundances of fecal ASV (A), genera (B), families (C), and phyla (D). The heritability value is represented by the dot and the SE by the line interval. The vertical dashed line illustrates the 0.10 heritability threshold. Only significant estimates are presented ($P \le 0.05$). Blue and red dots correspond to heritability estimates that passed and did not pass the FDR 10%, respectively.</figcaption></figure>
we found multiple significant genetic correlations: 63 between ASV (P between $6.38 \times 10^{-64}$ and 0.049), 58 between genera (P between $1.24 \times 10^{-93}$ and 0.050), 31 between families (P between $1.86 \times 10^{-12}$ and 0.049), and 7 between phyla (P between $1.43 \times 10^{-5}$ and 0.015). To account for multiple comparisons, we also applied a correction to the genetic correlations for the number of independent tests (i.e., 815 for ASV, 774 for genus, 186 for family, and 10 for phylum). Considering an FDR 10% threshold, 17 correlations between ASV remained significant, ranging in values from 0.682 (SE = 0.194) to values close to 1.0 (SE between 0.057 and 0.292). At the genus level, 20 genetic correlations were still significant. Seven of these correlations were negative (r between -0.99 ± 0.23 and -0.77 ± 0.21) and the 10 positive genetic correlations ranged from 0.80 (SE = 0.26) to estimates close to 1.0 (SE between 0.05 and 0.21). Specifically, the most significant correlation was between *Paeniclostridium* and *Romboutsia* (r close to 1.0, SE = 0.05, $1.24 \times 10^{-93}$, q = 9.59 × 10⁻⁹). At the upper taxonomic rank, 15 correlations between families remained. Four of these estimates were negative (r between -0.90 ± 0.31 and -0.65 ± 0.20), and the 11 others were positive with values ranging between 0.83 (SE = 0.17) and estimates close to 1.0 (SE between 0.14 and 0.28). At the highest level we studied, the 7 significant genetic correlations between phyla (P ≤ 0.50) still presented significant FDR 10% corrected P-values (q between $1.4 \times 10^{-4}$ and 0.21). Two of these genetic correlations were negative and both involved *Cyanobacteria* with *Bacteroidota* (r close to 1.0, SE = 0.24) and with *Patescibacteria* (r = -0.99, SE = 0.30). All other significant correlations ranged from 0.60 and value close to 1.0 (SE between 0.240 and 0.356).

Genetic correlations were also estimated between α- and β-diversity indices. Only NMDS1_Jaccard presented significant genetic correlations with NMDS1_Bray-Curtis (r = 0.852, SE = 0.135, P = $1.72 \times 10^{10}$) and Observed Richness (r = 0.852, SE = 0.135, P = $1.72 \times 10^{10}$). However, after correcting for multiple-testing, only the genetic correlation between the NMDS1_Jaccard and NMDS1_Bray-Curtis remained significant (q = $2.59 \times 10^{-9}$).

Furthermore, both α- and β-diversity indices presented significant phenotypic and genetic correlations with taxa abundances (Figure 4; Supplemental Figure S2; see Notes). Specifically, α-diversity indices were significantly genetically correlated with 3 ASV, 2 families, and one phylum. However, only the genetic correlation with the *Patescibacteria* phylum remained significant (r = -0.94, SE = 0.31) after applying the FDR 10% correction. Regarding β-diversity metrics, we observed 19 genetic correlations with ASV, 17 with genera, 15 with families,

and 7 with phyla. Interestingly, most genetic correlations between taxa abundances and β-diversity indices were negative: 16, 11, 14, and 6 negative genetic correlations were found for ASV, genus, family, and phylum taxonomic levels, respectively. After FDR 10% correction, 9 genetic correlations with ASV, 6 with genera, 10 with families, and 8 with phyla remained significant. Among these, only one significant estimate revealed a genetic association between the Shannon α-diversity index and the *Patescibacteria* phylum.

## Genetic Correlations Between Microbiota-Based Phenotypes and Cow Performance

In our analyses, we estimated genetic correlations between taxa abundances and host phenotypes (Figure 5). Eight ASV, 4 genera, 3 families, and 1 phylum showed significant genetic correlations (P ≤ 0.05) with at least 1 phenotype. The 8 ASV were involved in 9 significant correlations, consisting of 3 genetic correlations with production traits (MY, PY, and FC), 3 with fertility phenotypes (CCR and IFF), and 3 with SCS. The 4 genera with significant genetic correlations included: an unclassified genus from the *Peptococcaceae* family with MY (r = 0.50, SE = 0.25, P = 0.048), *Succinivibrio* with SCS (r = -0.85, SE = 0.38, P = 0.026), 2 unclassified genera from *Gastranaerophilales* and *WCHB1-41* orders with CCR (r = -0.94, SE = 0.30, P = 0.0015) and IFF (r = 0.92, SE = 0.42, P = 0.027), respectively. At the family level, 2 genetic correlations were found between SCS and *Succinivibrionaceae* (r = -0.45, SE = 0.21, P = 0.032) and *Erysipelatoclostridiaceae* (r = 0.46, SE = 0.23, P = 0.048). *Erysipelotrichaceae* was also genetically correlated with IFF (r = -0.91, SE = 0.34, P = 0.007) and *Peptococcaceae* was genetically correlated with MY (r = 0.53, SE = 0.25, P = 0.037). The only phylum to be genetically correlated with a host trait was *Proteobacteria* with SCS (r = -0.36, SE = 0.17, P = 0.038). Notably, SCS was the host phenotype most genetically correlated with taxa abundances. Nevertheless, after applying the FDR 10% correction, no genetic correlations remained significant between host traits and taxa abundances. Similar trends were observed for genetic correlations between host traits and diversity indices. However, IFF was genetically correlated with both NMDS2_Bray-Curtis (r = 0.88, SE = 0.41, P = 0.031) and NMDS2_Jaccard (r = 0.99, SE = 0.51, P = 0.049).

Conversely, although the estimated values of the phenotypic correlations were lower, a higher number of them presented significant correlations with taxa abundances and diversity indexes (Supplemental Figures S3 and S4; see Notes).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Genetic (upper) and phenotypic (lower) correlations between CLR-transformed abundances of heritable ASV (A), genus (B), family (C), and phylum (D) taxonomic ranks. The value of the correlations is represented by a color scale from blue (r = -1) to red (r = +1). Corr = correlation.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Genetic correlations between diversity indices and abundances of heritable ASV (A), genus (B), family (C), and phylum (D) taxonomic ranks. Diversity indices were calculated on a data set rarefied to 5,000 reads. Shannon and Observed Richness refer to α-diversity. NMDS1\_Bray-Curtis and NMDS2\_Bray-Curtis illustrated the values of the coordinates of the 2 first axes of NMDS constructed on the Bray-Curtis dissimilarity matrix. NMDS1\_Jaccard and NMDS2\_Jaccard illustrate the values of the coordinates of the 2 first axes of NMDS constructed on the Jaccard dissimilarity matrix. Taxa abundances referred to the CLR-transformed abundances of taxa. The value of the correlations is represented by a color scale from blue (r = −1) to red (r = +1). Corr = correlation.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Significant genetic correlations between host's phenotypes and CLR-transformed taxa abundances of heritable taxa ($h^2$ ± SE >0). Host's phenotypes referred to the YD of the host performance (MY = milk yield, PY = protein yield, FC = fat content, IFF = interval from first AI to fertilizing AI, CCR = conception rate) and are represented by yellow nodes. Taxa abundances referred to the CLR-transformed abundances of taxa. Taxa are represented with green nodes and the taxonomic rank is specified with dark green for ASV (lowest) and light green for phylum (highest) taxonomic ranks. The genetic correlations are represented with traits from blue (-1.0) to red (1.0).</figcaption></figure>

## DISCUSSION

Our investigation, based on a cohort of 1,875 Holstein cows sourced from 144 commercial herds, marks a major effort in deciphering the genetic parameters governing the fecal microbiota within milk production farms. Through our comprehensive approach, we successfully revealed the influence of host genetics on both the diversity and composition of the fecal microbiota, as well as relevant genetic correlations between microbiota characteristics and key traits of interest in Holstein cows.

### Insights into the Characterization of Fecal Microbiota in Dairy Cows

Our analysis reveals a rich and diverse fecal microbiota in lactating Holstein cows. In line with other studies, *Firmicutes* predominated, accounting for over half of the sequences (Mao et al., 2015; Cendron et al., 2020). At the genus level, our study identified an unclassified genus from the *UCG-010* family (*Oscillospirales* order) as the most abundant, contrasting with the findings of Cendron et al. (2020), who reported *Rikenellaceae RC9 gut group* and *Bacteroides* as the most abundant genera in fecal samples of Holstein cows. Despite discrepancies at the genus level, our study aligns with previous research on ruminant gastrointestinal microbiota, where the *Rikenellaceae RC9 gut group*'s associated family, *Rikenellaceae*, emerged as the most abundant family in the feces (Bach et al., 2019; Zhang et al., 2019; Cendron et al., 2020).

Interestingly, our study revealed many taxa with high prevalence, supporting the existence of a core microbiota shared among the majority of Holstein cows under study. This core microbiota observed in at least 90% of the sample, comprised 0.09% of the ASV, 11.2% of the genera, 17.8% of the families, and 24.2% of the phyla. Five genera, 7 families, and 3 phyla were consistently identified in all samples. At the phyla level, these findings are consistent with previous research reporting the prevalence of *Firmicutes*, *Bacteroidetes*, and *Proteobacteria* within the core microbiota of bovine gastrointestinal tracts (Mao et al., 2015; Gu et al., 2023). Additionally, our study corroborates previous observations with high prevalence of *Spirochaetota* and *Cyanobacteria* phyla within the core microbiota (Mao et al., 2015).

### Role of Host Genetics in Shaping the Fecal Microbiota Composition and Diversity

Our results confirm the hypothesis that the host genetics exert a role in the composition and diversity of the fecal microbiota in adult cows. Heritability estimates associated with taxa of the core microbiota suggest that the overall composition of the fecal microbiota in cows can be modulated by the host genetics.

Our analysis revealed low to moderate heritabilities for a significant proportion of the analyzed taxa, that is, 9% of the ASV, 16% of the genera, 21% of the families, and 27% of the phyla tested, with estimates ranging from 0.08 to 0.31. To our knowledge, similar results regarding
the fecal microbiota of cows have not been previously reported. However, despite differences between the ruminal and fecal microbiotas in ruminants (Mao et al., 2015; Monteiro et al., 2022), similar results were obtained in both compartments for which a notable proportion of nonheritable taxa is generally observed, with the remaining taxa showing low-to-moderate heritabilities. Indeed, in Holstein cows, various studies observed a large proportion (90%-94%) of nonheritable ruminal taxa, while some taxa exhibited heritabilities ranging from 0.08 to 0.48 for operational taxonomic units (OTU) and 0.17 to 0.41 for genera (Difford et al., 2018; Saborío-Montero et al., 2020; Zhang et al., 2020). In beef cattle, heritability estimates vary among studies, with values ranging from 0 to 0.25 for 56 bacterial taxa in the rumen (Li et al., 2019) to considerably higher estimates reaching 0.82 for OTU, families and phyla (Abbas et al., 2020). Similarly, in sheep, heritabilities of rumen bacterial abundances ranged from 0 to 0.29 (Martinez Boggio et al., 2022). Our findings are also consistent with observations of gut and fecal microbiota across species, with heritability estimates for taxa abundances ranging from 0.067 to 0.21 in human and baboon gut microbiota (Grieneisen et al., 2021) and a broad spectrum of 0.025 to 0.55 in pigs (Bergamaschi et al., 2020).

In our study, the most abundant taxa were not systematically associated with host genetics. For instance, the most abundant genus, *Rikenellaceae RC9 gut group*, was not found to be heritable. Similarly, although *Firmicutes* was the most abundant phylum, it was not found to be heritable in our analysis. This statement contradicts the observations made in the rumen of beef cattle, where *Firmicutes*, although less prevalent, was moderately heritable (h² = 0.30; Abbas et al., 2020). In contrast, *Proteobacteria* and *Bacteroidetes*, 2 other phyla present in all samples but at lower abundance than *Firmicutes*, were found to be the most heritable phyla, with estimates of 0.24 and 0.18, respectively. This aligns with findings in beef cattle rumen microbiota, where the most dominant phylum, *Bacteroides*, presented similar heritability values (h² = 0.23), and *Proteobacteria* was identified as the most heritable microbial feature (h² = 0.53; Abbas et al., 2020).

Our investigation into the heritability of microbial β-diversity, while revealing a moderate value (h² = 0.14), aligns with previous studies conducted on the rumen microbiota (from 0 to 0.28; Difford et al., 2018; Li et al., 2019). This finding emphasizes the influence of host genetics on the overall composition of the fecal microbiota. In contrast, our study did not find significant heritability for α-diversity indices, which is inconsistent with previous studies on rumen (Li et al., 2019) and monogastric species, including pigs (Déru et al., 2022) and humans (Grieneisen et al., 2021).

Echoing the observations of Difford et al. (2018), who also corrected for multiple tests in their rumen microbiota study, we found that only a limited number of taxa remained heritable after adjusting for multiple tests using the FDR 10% threshold. It is important to note that this multiple-testing correction method is known to be conservative, particularly in microbiota studies where the number of phenotyped individuals is often limited. Nevertheless, these preliminary findings in dairy cattle feces underscore the complexity of host-microbiota interactions. These interactions are largely influenced by environmental factors, such as the animal's diet, and the analytical methods employed. Despite these complexities, our study holds promise. Indeed, we were able to highlight significant low-to-moderate heritabilities for taxa abundances in the feces of Holstein lactating cows reared in commercial farms, even with unknown information regarding their diet, farm management, and other environmental management factors. Thus, these findings suggest that it may be possible to consider some taxa abundances as potential breeding objectives with large-scale analyses and commercial samplings.

## Genetic Correlations Between Fecal Microbiota Composition and Diversity

Genetic correlations across all taxonomic levels surpassed phenotypic correlations, aligning with established expectations (Cheverud, 1988). Various host environmental factors, such as parity or herd management practices, could potentially reduce phenotypic correlation values (Rothschild et al., 2018). The high genetic correlations observed, with some nearing -1.0 or 1.0, underscore the complex interplay and intricate network within the feces, in line with the findings of Zhang et al. (2019). Despite some imprecisions due to high SE in some genetic correlations, many passed the multiple-testing correction, indicating strong trends.

Interestingly, the phylogenetic proximity did not necessarily correlate with stronger genetic correlations. As an example, the most significant genetic correlations between taxa allowed for the identification of a notable cluster composed of *Paenicrostridium* (*Peptostreptococcaceae* family), *Romboutsia* (*Peptostreptococcaceae* family), *Turicibacter* (*Erysipelotrichaceae* family), and *Clostridium sensu stricto* 1 (*Clostridiaceae* family). Despite belonging to different families, these genera exhibited remarkably high genetic correlations ranging from 0.90 to 1.0. This result aligns with previous studies that have demonstrated correlations or co-occurrence of these genera in the ruminant gut (Wang et al., 2018; Kwon et al., 2021). Additionally, Castro-Ramos et al. (2022) have shown a close relationship between these genera during the anaerobic digestion of bovine manure,
all of which have been correlated with various volatile
fatty acids.

Our research unveils the complex networks and in-
teractions within the fecal community, suggesting that
selection on specific taxa may have broader implications
for the overall structure and composition of the microbial
community. This is exemplified by strong genetic correla-
tions between taxa in upper levels, such as *Bacteroidetes*
and *Proteobacteria* (r = 0.44), which ranked as the second
and fourth most abundant phyla in the feces of Holstein
cattle. This finding is further underscored by genetic
correlations between taxa and both α- and β-diversity
indices, implying that some of the most prevalent taxa
may play key roles in shaping the overall diversity and
structure of the fecal microbiota. This highlights their
importance in maintaining microbial community stability
and function.

Genetic Correlations Between the Fecal Microbiota
and Host Performances

Although some studies have explored genetic correla-
tions between the ruminal microbiota and certain traits of
interest in dairy cattle (Xue et al., 2018; Li et al., 2019;
Wallace et al., 2019; Gebreyesus et al., 2020; Saborío-
Montero et al., 2020), to the best of our knowledge, no
study has investigated genetic correlations between fecal
microbiota features and key phenotypes in this species.
Yet, stronger phenotypic associations were found between
milk production in dairy cows and their fecal microbiota
compared with their ruminal microbiota (Monteiro et
al., 2022). This underscores the importance of the fecal
microbiota in shaping key phenotypes in dairy cattle and
potentially enhancing the host's ability to metabolize a
larger proportion of its diet.

In the fecal microbiota, we identified genetic links
between certain taxa and dairy production. Notably, we
observed positive correlations between the *Peptococcaceae*
family and MY. This association may be attributed
to the butyrate production activity of *Peptococcaceae*,
known to enhance dairy production in cows (Urrutia et
al., 2019) and improve host's health, particularly through
anti-inflammatory properties (Fan et al., 2021). Ad-
ditionally, at a higher taxonomic level, we highlighted
various significant genetic correlations between host
traits and some taxa associated with the *Oscillospirales*
order. Among them, one ASV was negatively correlated
with both MY and PY, and 3 others demonstrated posi-
tive genetic correlations with FC, CCR, and SCS. These
results emphasize the importance of working at lower
taxonomic level to be able to detect the actions of spe-
cific taxa that could defer within the same taxonomic
group. Moreover, the examples of the *Peptococcaceae*
family and *Oscillospirales* order underscore the potential

role of some taxonomic groups in the expression of the
host's phenotype. It would be interesting to better under-
stand these genetic relationships, especially to find the
best strategy to implement microbiota data into predic-
tion equation, as proposed by Legarra and Christensen
(2023) or Martínez-Álvaro et al. (2022b). The use of a
specific taxonomic rank could be impactful, as suggested
by Zhang et al. (2020), but considering other methods
of microbiota analysis, such as the whole metagenomic
shotgun, could be interesting for working at the level
of microbial genes, as done by Martínez-Álvaro et al.
(2022b).

Some host nonproduction related phenotypes, such as methane emission, have been found to be strongly linked with cattle microbiota with genetic correlations varying between -0.68 and 0.69 (Difford et al., 2018; Saborío-Montero et al., 2020; Zhang et al., 2020; Martínez-Álvaro et al., 2022a). Our results also revealed genetic correlations between fecal taxa and udder health, as assessed by SCS, aligning with evidence suggesting the active role of the ruminant feces in immune functions (O'Hara et al., 2020). Previous studies have highlighted associations between the gut microbiota and host health, particularly in calves (Du et al., 2023), but such associations have been less explored in adults. We estimated negative genetic correlations between SCS and genera such as *Succinivibrio*, as well as its respective family and phylum, *Succinivibrionaceae* and *Proteobacteria*. These taxa were negatively correlated with SCS, as previously suggested (Zhong et al., 2018). Additionally, SCS was genetically correlated with 2 ASV, another genus, and another family. Thus, SCS and other health traits should also be considered as phenotypes that can be partly predicted with microbiota data. Aside from health-related traits, fertility traits also presented some significant genetic correlations with taxa abundances, especially with 3 ASV, one genus, and one phylum. Interestingly, one of the fertility traits, IFF, was also the only host's YD to be correlated with β-diversity indices (i.e., NMDS2_Bray-Curtis and NMDS2_Jaccard). These findings suggested a potential interaction between gut microbiota and fertility, a phenomenon also observed in humans (Patel et al., 2022).

The use of YD in this study, which are performances
adjusted with accurate nongenetic effects estimated at
the national level, may explain the reduced correlations
observed, because the phenotypes compared the whole
cows' career with a sampling defining the fecal micro-
biota at a certain time. Furthermore, none of these low
genetic correlations between microbiota-based pheno-
types and host traits passed the rigorous multiple-testing
correlations, partly due to high SE. Despite this, these
findings provide valuable insights into potential genetic
associations between fecal microbiota composition and
host performance. Considering these associations could
be beneficial in a breeding context, because selecting for a healthier fecal microbiota could indirectly influence host performance (Zhang et al., 2020; Martínez-Álvaro et al., 2022b). Moreover, they highlight the potential precision gain achievable by incorporating fecal microbiota data into prediction models, as suggested by Legarra and Christensen (2023).

## CONCLUSIONS

This investigation of the fecal microbiota of Holstein cows from commercial herds provides insights into its structure and the influence of host genetics, with low to moderate heritability estimates of taxa. Our analyses reveal moderate to high genetic correlations between microbiota taxa, highlighting the complex network in Holstein feces. Future research should further explore genomic associations with microbiota. Identifying genes and metabolic pathways involved in the genetic determination could unravel the interplay between host genetics and microbial communities, enhancing our understanding of host-microbiotic interactions in dairy cattle. Such analyses may allow for the consideration of fecal microbiota as a breeding goal. Genetic correlations between abundances of fecal taxa and key phenotypes suggest promising applications in selection. More detailed analyses could help determine the extent to which prediction accuracy could be improved by incorporating microbiota information into selection equations. Fecal microbiota could then be considered as a large-scale phenotype and data source for improving animal health and optimizing production efficiency.

## NOTES

This work has been financially supported by the French "Association Nationale de la Recherche et de la Technologie" (ANRT, Paris, France) that allows the PhD contract of Louise Brulin (CIFRE n°2021/0713). The authors would like to thank the breeders and the technicians who were involved in the collection of the samples, as well as Gérard Heumez (Gènes Diffusion, Douai, France) for his assistance in the data processing. The authors gratefully acknowledge the French National Database (CTIG, Jouy-en-Josas, France) and Valogène (Paris, France) for providing phenotype and genotype data, respectively; and the contribution of the 1,000 Bull Genomes Project. The authors also thank the INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020. MIGALE bioinformatics facility, Jouy-en-Josas, France; https://doi.org/10.15454/1.5572390655343293E12) for providing the computational resources needed for our analyses. GD Biotech expresses gratitude to the ANR-10-EQPX-07-01 LIGAN-PM-UMR8199 (Lille, France)

for its technical support regarding the MiSeq sequencing. Supplemental material for this article is available at https://doi.org/10.6084/m9.fishshare.25534660. This work was conducted on farm animals reared for commercial purposes in compliance with the French regulation (Code Rural et de la Pêche Maritime) and the European Council Directive 98/58/EC. In accordance with the legislation, no approval from the Institutional Animal Care and Use Committee or ethics committee was necessary, as all the performed sampling operations were part of routine animal manipulations by duly authorized technicians. These technicians were employees of the Gènes Diffusion breeding company (Douai, France), who took part in the routine reproduction monitoring. They were holders of the CAFTI diploma (Certificat d'Aptitude aux Fonctions de Technicien d'Insémination – Certificate of Fitness for Insemination Technician Functions), approved by the declaration to the EDE (Departmental Establishments of Breeding), which authorizes them for biological sampling in agreement with animal welfare regulation. The farmers involved in this study agreed to the use of their animals' samples for research purposes. Authors LB, SD, GE, SMa, SMe, and CA were employed by GD Biotech/Gènes Diffusion company. The authors have not stated any other conflicts of interest.

**Nonstandard abbreviations used:** ASV = amplicon sequence variant; CAI<sub>f</sub> = interval between calving and first AI; CCR = conception rate; CLR = centered-log ratio; FC = fat content; FDR = false discovery rate; FDR 10% = FDR threshold of 10%; FY = fat yield; IFF = interval between first AI and fertilizing AI; MY = milk yield; NMDS = nonmetric multidimensional scaling; OTU = operational taxonomic unit; PC = protein content; PY = protein yield; RA = relative abundance; YD = yield deviation.

## REFERENCES

Abbas, W., J. T. Howard, H. A. Paz, K. E. Hales, J. E. Wells, L. A. Kuehn, G. E. Erickson, M. L. Spangler, and S. C. Fernando. 2020. Influence of host genetics in shaping the rumen bacterial community in beef cattle. Sci. Rep. 10:15101. https://doi.org/10.1038/s41598-020-72011-9.

Andrews, S. 2010. FastQC: A Quality Control Tool for High Throughput Sequence Data. Babraham Bioinformatics. Accessed Dec. 22, 2022. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/.

Bach, A., A. López-García, O. González-Recio, G. Elcoso, F. Fábregas, F. Chaucheyras-Durand, and M. Castex. 2019. Changes in the rumen and colon microbiota and effects of live yeast dietary supplementation during the transition from the dry period to lactation of dairy cows. J. Dairy Sci. 102:6180-6198. https://doi.org/10.3168/jds.2018-16105.

Bergamaschi, M., C. Maltecca, C. Schillebeeckx, N. P. McNulty, C. Schwab, C. Shull, J. Fix, and F. Tiezzi. 2020. Heritability and genome-wide association of swine gut microbiome features with growth and fatness parameters. Sci. Rep. 10:10134. https://doi.org/10.1038/s41598-020-66791-3.
Bozzi, D., J. A. Rasmussen, C. Carøe, H. Sveier, K. Nordøy, M. T. P. Gilbert, and M. T. Limborg. 2021. Salmon gut microbiota correlates with disease infection status: potential for monitoring health in farmed animals. Anim. Microbiome 3:30. https://doi.org/10.1186/s42523-021-00096-2.

Callahan, B. J., P. J. McMurdie, M. J. Rosen, A. W. Han, A. J. A. Johnson, and S. P. Holmes. 2016. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13:581–583. https://doi.org/10.1038/nmeth.3869.

Castro-Ramos, J. J., A. Solís-Oba, M. Solís-Oba, C. L. Calderón-Vázquez, J. M. Higuera-Rubio, and R. Castro-Rivera. 2022. Effect of the initial pH on the anaerobic digestion process of dairy cattle manure. AMB Express 12:162. https://doi.org/10.1186/s13568-022-01486-8.

Cendron, F., G. Niero, G. Carlino, M. Penasa, and M. Cassandro. 2020. Characterizing the fecal bacteria and archaea community of heifers and lactating cows through 16S rRNA next-generation sequencing. J. Appl. Genet. 61:593–605. https://doi.org/10.1007/s13353-020-00575-3.

Chen, Y., M. Akhtar, Z. Ma, T. Hu, Q. Liu, H. Pan, X. Zhang, A. A. Nafady, A. R. Ansari, E.-S. M. Abdel-Kafy, D. Shi, and H. Liu. 2023. Chicken cecal microbiota reduces abdominal fat deposition by regulating fat metabolism. NPJ Biofilms Microbiomes 9:1–16. https://doi.org/10.1038/s41522-023-00390-8.

Cheverud, J. M. 1988. A comparison of genetic and phenotypic correlations. Evolution 42:958–968. https://doi.org/10.2307/2408911.

Déru, V., A. Bouquet, O. Zemb, B. Blanchet, M. L. De Almeida, L. Cauquil, C. Carillier-Jacquin, and H. Gilbert. 2022. Genetic relationships between efficiency traits and gut microbiota traits in growing pigs being fed with a conventional or a high-fiber diet. J. Anim. Sci. 100:skac183. https://doi.org/10.1093/jas/skac183.

Difford, G. F., D. R. Plichta, P. Løvendahl, J. Lassen, S. J. Noel, O. Højberg, A.-D. G. Wright, Z. Zhu, L. Kristensen, H. B. Nielsen, B. Guldbrandtsen, and G. Sahana. 2018. Host genetics and the rumen microbiome jointly associate with methane emissions in dairy cows. PLoS Genet. 14:e1007580. https://doi.org/10.1371/journal.pgen.1007580.

Du, Y., Y. Gao, M. Hu, J. Hou, L. Yang, X. Wang, W. Du, J. Liu, and Q. Xu. 2023. Colonization and development of the gut microbiome in calves. J. Anim. Sci. Biotechnol. 14:46. https://doi.org/10.1186/s40104-023-00856-x.

Fan, P., C. D. Nelson, J. D. Driver, M. A. Elzo, F. Peñagaricano, and K. C. Jeong. 2021. Host genetics exerts lifelong effects upon hindgut microbiota and its association with bovine growth and immunity. ISME J. 15:2306–2321. https://doi.org/10.1038/s41396-021-00925-x.

Gebreyesus, G., G. F. Difford, B. Buitenhuis, J. Lassen, S. J. Noel, O. Højberg, D. R. Plichta, Z. Zhu, N. A. Poulsen, U. K. Sundekilde, P. Lovendahl, and G. Sahana. 2020. Predictive ability of host genetics and rumen microbiome for subclinical ketosis. J. Dairy Sci. 103:4557–4569. https://doi.org/10.3168/jds.2019-17824.

Grieneisen, L., M. Dasari, T. J. Gould, J. R. Björk, J.-C. Grenier, V. Yotova, D. Jansen, N. Gottel, J. B. Gordon, N. H. Learn, L. R. Gesquiere, T. L. Wango, R. S. Mututua, J. K. Warutere, L. Siodi, J. A. Gilbert, L. B. Barreiro, S. C. Alberts, J. Tung, E. A. Archie, and R. Blekhman. 2021. Gut microbiome heritability is nearly universal but environmentally contingent. Science 373:181–186. https://doi.org/10.1126/science.aba5483.

Gu, F., S. Zhu, J. Hou, Y. Tang, J.-X. Liu, Q. Xu, and H.-Z. Sun. 2023. The hindgut microbiome contributes to host oxidative stress in postpartum dairy cows by affecting glutathione synthesis process. Microbiome 11:87. https://doi.org/10.1186/s40168-023-01535-9.

Hill, J., C. McSweeney, A.-D. G. Wright, G. Bishop-Hurley, and K. Kalantar-zadeh. 2016. Measuring methane production from ruminants. Trends Biotechnol. 34:26–35. https://doi.org/10.1016/j.tibtech.2015.10.004.

Jami, E., B. A. White, and I. Mizrahi. 2014. Potential role of the bovine rumen microbiome in modulating milk composition and feed efficiency. PLoS One 9:e85423. https://doi.org/10.1371/journal.pone.0085423.

Kwon, M.-S., H. E. Jo, J. Lee, K.-S. Choi, D. Yu, Y. Oh, J. Park, and H.-J. Choi. 2021. Alteration of the gut microbiota in post-weaned calves following recovery from bovine coronavirus-mediated diarrhea. J. Anim. Sci. Technol. 63:125–136. https://doi.org/10.5187/jast.2021.e20.

Lahti, L., and S. Shetty. 2017. Tools for microbiome analysis in R. Version 1.24.0. Accessed Sep. 23, 2024. https://github.com/microbiome/microbiome.

Legarra, A., and O. F. Christensen. 2023. Genomic evaluation methods to include intermediate correlated features such as high-throughput or omics phenotypes. JDS Commun. 4:55–60. https://doi.org/10.3168/jdsc.2022-0276.

Li, F., C. Li, Y. Chen, J. Liu, C. Zhang, B. Irving, C. Fitzsimmons, G. Plastow, and L. L. Guan. 2019. Host genetics influence the rumen microbiota and heritable rumen microbial features associate with feed efficiency in cattle. Microbiome 7:92. https://doi.org/10.1186/s40168-019-0699-1.

Mao, S., M. Zhang, J. Liu, and W. Zhu. 2015. Characterising the bacterial microbiota across the gastrointestinal tracts of dairy cattle: Membership and potential function. Sci. Rep. 5:16116. https://doi.org/10.1038/srep16116.

Martínez-Álvaro, M., M. D. Auffret, C.-A. Duthie, R. J. Dewhurst, M. A. Cleveland, M. Watson, and R. Roehe. 2022a. Bovine host genome acts on rumen microbiome function linked to methane emissions. Commun. Biol. 5:350. https://doi.org/10.1038/s42003-022-03293-0.

Martínez-Álvaro, M., J. Mattock, M. Auffret, Z. Weng, C.-A. Duthie, R. J. Dewhurst, M. A. Cleveland, M. Watson, and R. Roehe. 2022b. Microbiome-driven breeding strategy potentially improves beef fatty acid profile benefiting human health and reduces methane emissions. Microbiome 10:166. https://doi.org/10.1186/s40168-022-01352-6.

Martinez Boggio, G., A. Meynadier, A. J. Buitenhuis, and C. Marie-Etancelin. 2022. Host genetic control on rumen microbiota and its impact on dairy traits in sheep. Genet. Sel. Evol. 54:77. https://doi.org/10.1186/s12711-022-00769-9.

McMurdie, P. J., and S. Holmes. 2013. phyloseq: An R Package for Reproducible Interactive Analysis and Graphics of Microbiome Census Data. PLoS One 8:e61217. https://doi.org/10.1371/journal.pone.0061217.

Meyer, K. 2007. WOMBAT—A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML). J. Zhejiang Univ. Sci. B 8:815–821. https://doi.org/10.1631/jzus.2007.B0815.

Monteiro, H. F., Z. Zhou, M. S. Gomes, P. M. G. Peixoto, E. C. R. Bonaglia, I. F. Canisso, B. C. Weimer, and F. S. Lima. 2022. Rumen and lower gut microbiomes relationship with feed efficiency and production traits throughout the lactation of Holstein dairy cows. Sci. Rep. 12:4904. https://doi.org/10.1038/s41598-022-08761-5.

Na, S. W., and L. L. Guan. 2022. Understanding the role of rumen epithelial host-microbe interactions in cattle feed efficiency. Anim. Nutr. 10:41–53. https://doi.org/10.1016/j.aninu.2022.04.002.

O'Hara, E., A. L. A. Neves, Y. Song, and L. L. Guan. 2020. The role of the gut microbiome in cattle production and health: Driver or passenger? Annu. Rev. Anim. Biosci. 8:199–220. https://doi.org/10.1146/annurev-animal-021419-083952.

Oksanen, J., G. L. Simpson, F. G. Blanchet, R. Kindt, P. Legendre, P. R. Minchin, R. B. O'Hara, P. Solymos, M. H. H. Stevens, E. Szoecs, H. Wagner, M. Barbour, M. Bedward, B. Bolker, D. Borcard, G. Carvalho, M. Chirico, M. D. Caceres, S. Durand, H. B. A. Evangelista, R. FitzJohn, M. Friendly, B. Furneaux, G. Hannigan, M. O. Hill, L. Lahti, D. McGlinn, M.-H. Ouellette, E. R. Cunha, T. Smith, A. Stier, C. J. F. T. Braak, and J. Weedon. 2022. vegan: Community Ecology Package.

Patel, N., N. Patel, S. Pal, N. Nathani, R. Pandit, M. Patel, N. Patel, C. Joshi, and B. Parekh. 2022. Distinct gut and vaginal microbiota profile in women with recurrent implantation failure and unexplained infertility. BMC Womens Health 22:113. https://doi.org/10.1186/s12905-022-01681-6.

Quast, C., E. Pruesse, P. Yilmaz, J. Gerken, T. Schweer, P. Yarza, J. Peplies, and F. O. Glöckner. 2013. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41:D590–D596. https://doi.org/10.1093/nar/gks1219.
R Core Team. 2022. R: A Language and Environment for Statistical Computing. Version 4.3.2. R Foundation for Statistical Computing, Vienna, Austria.

Reyer, H., M. Oster, U. M. McCormack, E. Muráni, G. E. Gardiner, S. Ponsuksili, P. G. Lawlor, and K. Wimmers. 2020. Host-microbiota interactions in ileum and caecum of pigs divergent in feed efficiency contribute to nutrient utilization. Microorganisms 8:563. https://doi.org/10.3390/microorganisms8040563.

Rothschild, D., O. Weissbrod, E. Barkan, A. Kurilshikov, T. Korem, D. Zeevi, P. I. Costea, A. Godneva, I. N. Kalka, N. Bar, S. Shilo, D. Lador, A. V. Vila, N. Zmora, M. Pevsner-Fischer, D. Israeli, N. Kosower, G. Malka, B. C. Wolf, T. Avnit-Sagi, M. Lotan-Pompan, A. Weinberger, Z. Halpern, S. Carmi, J. Fu, C. Wijmenga, A. Zhernakova, E. Elinav, and E. Segal. 2018. Environment dominates over host genetics in shaping human gut microbiota. Nature 555:210-215. https://doi.org/10.1038/nature25973.

Saborío-Montero, A., M. Gutiérrez-Rivas, A. García-Rodríguez, R. Atxaerandio, I. Goiri, E. López de Maturana, J. A. Jiménez-Montero, R. Alenda, and O. González-Recio. 2020. Structural equation models to disentangle the biological relationship between microbiota and complex traits: Methane production in dairy cattle as a case of study. J. Anim. Breed. Genet. 137:36-48. https://doi.org/10.1111/jbg.12444.

Urrutia, N., R. Bomberger, C. Matamoros, and K. J. Harvatine. 2019. Effect of dietary supplementation of sodium acetate and calcium butyrate on milk fat synthesis in lactating dairy cows. J. Dairy Sci. 102:5172-5181. https://doi.org/10.3168/jds.2018-16024.

Wallace, R. J., G. Sasson, P. C. Garnsworthy, I. Tapio, E. Gregson, P. Bani, P. Huhtanen, A. R. Bayat, F. Strozzi, F. Biscarini, T. J. Snelling, N. Saunders, S. L. Potterton, J. Craigon, A. Minuti, E. Trevisi, M. L. Callegari, F. P. Cappelli, E. H. Cabezas-Garcia, J. Vilkki, C. Pinares-Patino, K. O. Fliegerová, J. Mrázek, H. Sechovcová, J. Kopečný, A. Bonin, F. Boyer, P. Taberlet, F. Kokou, E. Halperin, J. L. Williams, K. J. Shingfield, and I. Mizrahi. 2019. A heritable subset of the core rumen microbiome dictates dairy cow productivity and emissions. Sci. Adv. 5:eaav8391. https://doi.org/10.1126/sciadv.aav8391.

Wang, Y., H. Zhang, L. Zhu, Y. Xu, N. Liu, X. Sun, L. Hu, H. Huang, K. Wei, and R. Zhu. 2018. Dynamic distribution of gut microbiota in

goats at different ages and health states. Front. Microbiol. 9:2509. https://doi.org/10.3389/fmicb.2018.02509.

Xu, Q., Q. Qiao, Y. Gao, J. Hou, M. Hu, Y. Du, K. Zhao, and X. Li. 2021. Gut microbiota and their role in health and metabolic disease of dairy cow. Front. Nutr. 8:701511. https://doi.org/10.3389/fnut.2021.701511.

Xue, M., H. Sun, X. Wu, L. L. Guan, and J. Liu. 2018. Assessment of rumen microbiota from a large dairy cattle cohort reveals the pan and core bacteriomes contributing to varied phenotypes. Appl. Environ. Microbiol. 84:e00970-18. https://doi.org/10.1128/AEM.00970-18.

Zhang, G., Y. Wang, H. Luo, W. Qiu, H. Zhang, L. Hu, Y. Wang, G. Dong, and G. Guo. 2019. The association between inflammaging and age-related changes in the ruminal and fecal microbiota among lactating Holstein cows. Front. Microbiol. 10:1803. https://doi.org/10.3389/fmicb.2019.01803.

Zhang, Q., G. Difford, G. Sahana, P. Løvendahl, J. Lassen, M. S. Lund, B. Guldbrandtsen, and L. Janss. 2020. Bayesian modeling reveals host genetics associated with rumen microbiota jointly influence methane emission in dairy cows. ISME J. 14:2019-2033. https://doi.org/10.1038/s41396-020-0663-x.

Zhong, Y., M. Xue, and J. Liu. 2018. Composition of rumen bacterial community in dairy cows with different levels of somatic cell counts. Front. Microbiol. 9:3217. https://doi.org/10.3389/fmicb.2018.03217.

Zhou, Q., F. Lan, S. Gu, G. Li, G. Wu, Y. Yan, X. Li, J. Jin, C. Wen, C. Sun, and N. Yang. 2023. Genetic and microbiome analysis of feed efficiency in laying hens. Poult. Sci. 102:102393. https://doi.org/10.1016/j.psj.2022.102393.

## ORCIDS

L. Brulin ◎ https://orcid.org/0009-0006-8264-1064

S. Ducrocq ◎ https://orcid.org/0009-0004-6898-156X

C. Audebert ◎ https://orcid.org/0000-0002-6509-7574

M.-P. Sanchez ◎ https://orcid.org/0000-0002-1371-5342