# Microbiability and microbiome-wide association analyses of feed efficiency and performance traits in pigs

Amir Aliakbari, Olivier Zemb, Laurent Cauquil, Céline Barilly, Yvon Billon, Hélène Gilbert

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

## To cite this version:

Amir Aliakbari, Olivier Zemb, Laurent Cauquil, Céline Barilly, Yvon Billon, et al.. Microbiability and microbiome-wide association analyses of feed efficiency and performance traits in pigs. Genetics Selection Evolution, 2022, 54 (1), 15 p. 10.1186/s12711-022-00717-7. hal-03650709

HAL Id: hal-03650709

https://hal.inrae.fr/hal-03650709

Submitted on 26 Apr 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Microbiability and microbiome-wide association analyses of feed efficiency and performance traits in pigs

Amir Aliakbari¹*, Olivier Zemb¹, Laurent Cauquil¹, Céline Barilly¹, Yvon Billon² and Hélène Gilbert¹

**Abstract**

**Background:** The objective of the present study was to investigate how variation in the faecal microbial composition is associated with variation in average daily gain (ADG), backfat thickness (BFT), daily feed intake (DFI), feed conversion ratio (FCR), and residual feed intake (RFI), using data from two experimental pig lines that were divergent for feed efficiency. Estimates of microbiability were obtained by a Bayesian approach using animal mixed models. Microbiome-wide association analyses (MWAS) were conducted by single-operational taxonomic units (OTU) regression and by back-solving solutions of best linear unbiased prediction using a microbiome covariance matrix. In addition, accuracy of microbiome predictions of phenotypes using the microbiome covariance matrix was evaluated.

**Results:** Estimates of heritability ranged from 0.31 ± 0.13 for FCR to 0.51 ± 0.10 for BFT. Estimates of microbiability were lower than those of heritability for all traits and were 0.11 ± 0.09 for RFI, 0.20 ± 0.11 for FCR, 0.04 ± 0.03 for DFI, 0.03 ± 0.03 for ADG, and 0.02 ± 0.03 for BFT. Bivariate analyses showed a high microbial correlation of 0.70 ± 0.34 between RFI and FCR. The two approaches used for MWAS showed similar results. Overall, eight OTU with significant or suggestive effects on the five traits were identified. They belonged to the genera and families that are mainly involved in producing short-chain fatty acids and digestive enzymes. Prediction accuracy of phenotypes using a full model including the genetic and microbiota components ranged from 0.60 ± 0.19 to 0.78 ± 0.05. Similar accuracies of predictions of the microbial component were observed using models that did or did not include an additive animal effect, suggesting no interaction with the genetic effect.

**Conclusions:** Our results showed substantial associations of the faecal microbiome with feed efficiency related traits but negligible effects with growth traits. Microbiome data incorporated as a covariance matrix can be used to predict phenotypes of animals that do not (yet) have phenotypic information. Connecting breeding environment between training sets and predicted populations could be necessary to obtain reliable microbiome predictions.

**Background**

Recent advances in the collection of microbiota information make it possible to analyze the interplay between complex traits and the microbial community of the gastrointestinal tract (GIT) in animals and humans. This

is especially essential in the pig industry since previous studies have suggested that the GIT microbiome could contribute to the variability of feed efficiency in pigs [1-3]. From the perspective of quantitative genetics, the effect of the microbiome on a trait can be quantified by estimating the microbiability [4], which is the proportion of phenotypic variance of the trait that is explained by between-animal differences in the microbial community. Estimation of the microbiability requires a microbial relationship matrix between host animals [4]. Using

<footnote>*Correspondence: amir.aliakbari@inrae.fr
¹GenPhySE, Université de Toulouse, INRAE, ENVT, 31320 Castanet-Tolosan, France
Full list of author information is available at the end of the article</footnote>
this approach, Camarinha-Silva et al. [3] reported higher
estimates of microbiability for feed conversion ratio
(FCR) (0.21±0.14) and feed intake (0.16±0.10) than
their corresponding heritabilities. Similarly, a recent
study revealed that the proportion of variance captured
by the microbiome for meat quality and carcass compo-
sition traits in crossbred pigs varied over time, with an
increased proportion from weaning to off-test for most
of the studied traits [5]. This study also reported micro-
biability (m²) estimates that were higher than the cor-
responding heritability (h²) estimates for some traits,
particularly at the off-test stage. In contrast, Tang et al.
[6] obtained a lower estimate of m² than that of h² for
body weight (BW), average daily gain (ADG), backfat
thickness (BFT), and intramuscular fatness, using sam-
ples from five different parts of the GIT. Overall, these
studies highlight the relevance of variability in the GIT
microbiome composition associated with variability in
performance traits, which suggests the possibility of pre-
dicting future phenotypes based on predicted microbial
values (m̂) [7]. However, in livestock, only a few studies
have evaluated the accuracy of phenotype predictions by
including microbiota effects in linear mixed models [7,
8]. In addition, similar to genome-wide association stud-
ies, microbiome components can be considered as poten-
tial markers of the selected complex traits, and their
associations can be identified through microbiome-wide
association studies (MWAS) [10]. In an early MWAS in
the Piétrain pig breed, Camarinha-Silva et al. [3] identi-
fied several operational taxonomic units (OTU) with a
relative abundance that was significantly associated with
ADG, FCR, or feed intake, and they suggested, by anal-
ogy to the polygenic determinism of traits, that these
traits could have a polymicrobial nature [3]. To the best
of our knowledge, there is no additional published litera-
ture on this topic in pigs, although there are numerous
examples in humans and a few in other livestock spe-
cies [10, 11]. Given these previous results, we hypoth-
esized that microbiota information can help to predict
and manage traits of interest in pigs. The first objective
of the present study was to investigate the association of
faecal microbial variants with feed efficiency (FCR and
residual feed intake (RFI)) and other performance traits,
including ADG, BFT, and daily feed intake (DFI). The
second objective was to evaluate the accuracy of micro-
biome predictions of phenotypes of animals that do not
(yet) have phenotypic information with the incorporation
of a microbiome covariance matrix. We used data from
two experimental pig lines that were divergently selected
for RFI, which ensured sizeable variability of the traits of
interest. In this work, we obtained microbiability esti-
mates for the traits using mixed models with or without
accounting for the genetic background of the hosts.

**Methods**

Population structure, studied traits and sampling

Phenotypic records and faecal samples (604 pigs) were
collected from the two last generations (G9 and G10) of
two experimental French Large White pig lines that were
developed over 11 generations of divergent selection for
RFI during 18 years at INRAE (UE GenESI, Surgères,
France, https://doi.org/10.15454/1.557241548118584
7E12), as described in Aliakbari et al. [12]. The struc-
ture of the lines and the selection process are described
in Aliakbari et al. [13]. After weaning (28 days of age),
all pigs were penned in groups of 24, by line and sex. At
10 weeks of age, pigs from each pen were distributed
in two growing-finishing pens (n=12 per pen). There
were four pens per room and one or two rooms per con-
temporary group (CG). Growing-finishing pens were
equipped with single-place electronic feeders ACEMA64
(ACEMO, France) to record individual feed intake. A pel-
leted diet based on cereals and soya bean meal was avail-
able ad libitum and contained 10MJ net energy (NE)/kg
and 160 g crude protein/kg, with a minimum of 0.80 g
digestible Lys/MJ NE. Animals had free access to water
at all stages. Complete pedigree information was regis-
tered, starting at least one generation before F0 ancestors
(founder generation of the divergent lines) until G10.

A set of 157 animals from generation G9 (entire male
candidates for selection) had records for feed intake
(DFI), feed efficiency traits (FCR and RFI), growth rate
(ADG) from 35 to 95 kg of body weight (BW), and live
body composition traits at 95 kg BW (BFT). The remain-
ing animals, from G10 (females and castrated males), had
records for growth rate, feed efficiency from 10 weeks of
age until slaughter (115 kg BW), and backfat thickness at
23 weeks of age (BFT23). Different multiple linear regres-
sion equations were used to compute realized RFI for
selection candidates versus G10 females and castrated
males, considering their test differences, as described
in [14]. First, the realized RFI for a selection candidate
was defined as the residual for that animal of a multi-
ple regression across selection candidates from G0 to
G9 of DFI on ADG and BFT (measured by ultrasound),
including the fixed effects of pen size and CG. Then, for
all females and castrated males from G1 to G10, the real-
ized RFI was the residual of multiple regression of DFI
on average metabolic body weight, ADG, and indica-
tors of body composition (carcass BFT and lean meat
content (computed from cut weights) until G9, BFT23
for G10 animals), including the fixed effects of pen size,
CG and sex and covariate of BW at the beginning of the
test. The FCR was computed based on the correspond-
ing test period of the two groups of animals. Standard-
ized phenotypes of RFI, FCR, DFI, ADG, and BFT were
used, as previously proposed in Aliakbari et al. [13]. The
descriptive information for these five traits and the sampled G9 and G10 animals is in Table 1.

The faecal samples were collected at 15 weeks of age to obtain the gut microbial information. Immediately after collection, the samples were homogenized and placed in dry ice, before storage at – 80 °C until DNA extraction (see next section).

## Microbial information

The Quick-DNA™ Faecal Microbe Miniprep Kit™ (Zymo Research, Freiburg, Germany) was used to extract microbial DNA based on a 15-min bead-beating step at 30 Hz. PCR amplification of the V3-V4 region of the 16S rRNA gene obtained from diluted genomic DNA was carried out by using two primers: F343 (CTTTCCCTACAC GACGCTCTTCCGATCTTACGGRAGGCAGCAG) and R784 (GGAGTTCAGACGTGTGCTCTTCCGATC TTACCAGGGTATCTAATCCT), in 30 cycles with an annealing temperature of 65 °C. To assemble pair-end sequences, the Flash software v1.2.6 [15] was used with an overlap of at least 10-bp between the forward and reverse sequences and by allowing 10% mismatches. Single multiplexing was performed using an in-house 6-bp index, which was added to R784 during a second PCR with 12 cycles and the forward primer (AATGATACG GCGACCACCGAGATCTACACTCTTTCCCTACAC GAC) and reverse primer (CAAGCAGAAGACGGC ATACGAGAT-index-GTGACTGGAGTTCAGACG TGT). Then, the resulting PCR products were purified and loaded on the Illumina MiSeq cartridge according to the manufacturer's instructions. The quality of the runs was checked internally using PhiX, and each pair-end sequence was assigned to its sample using the integrated index, with the bcl2fastq Illumina software. The sequences were submitted to the Bioproject database (https://www.ncbi.nlm.nih.gov/bioproject/) with accession number PRJNA701065. Reads were filtered and trimmed for high-quality sequences using the DADA2 package in the R software [16] with the following parameters: `maxN=0`, `maxEE=2`, `truncQ=2`, `trimleft=17`.

<table><caption>**Table 1** Descriptive statistics of the standardized phenotypes for the evaluated production traits</caption><thead><tr><th>Trait</th><th>Number</th><th>Min</th><th>Max</th><th>Average</th><th>SD</th></tr></thead><tbody><tr><td>RFI</td><td>522</td><td>– 0.38</td><td>0.39</td><td>0.00</td><td>0.15</td></tr><tr><td>FCR</td><td>548</td><td>1.60</td><td>3.93</td><td>2.78</td><td>0.33</td></tr><tr><td>DFI</td><td>542</td><td>1.37</td><td>2.95</td><td>2.20</td><td>0.29</td></tr><tr><td>ADG</td><td>575</td><td>0.51</td><td>1.01</td><td>0.78</td><td>0.08</td></tr><tr><td>BFT</td><td>541</td><td>9.8</td><td>46.6</td><td>23.3</td><td>10.0</td></tr></tbody></table>

RFI residual feed intake (kg/day), FCR feed conversion ratio (kg/kg), DFI daily feed intake (kg/day), ADG average daily gain (kg/day), BFT backfat thickness (mm)

Chimera were removed with the consensus option in DADA2 to obtain the final OTU abundance table. Since no further clustering was applied, the OTU were equivalent to amplicon sequence variants (ASV) in this study. Then, taxonomic annotation was predicted using the assignTaxonomy function of the DADA2 package with the Silva Dataset v132 [17].

To account for differences in sequencing depth across samples, a rarefaction step was applied: for all samples with more than 9000 reads, an equal number of 9000 reads were randomly selected. Sixteen samples had a smaller number of reads and were discarded for later analyses. After this rarefaction step, the final abundance table contained 5689 OTU for 588 samples (295 LRFI and 293 HRFI). Finally, following Rothschild et al. [18], OTU in the rarefied table were filtered for more than 1% non-zero values across all sampled animals, which decreased the number of OTU to 2630. Then, a constant value of 1 was added to the rarefied table, which allowed calculation of log values for OTU with zero abundances in the downstream analyses.

## Statistical analyses

### Estimation of variance components

For each trait, the following four univariate linear models were fitted:

$$y = Xb + e, \quad (1)$$

$$y = Xb + Z_2m + e, \quad (2)$$

$$y = Xb + Z_1a + e, \quad (3)$$

$$y = Xb + Z_1a + Z_2m + e, \quad (4)$$

where *y* is the vector of observations of each of the five traits, *b* is the vector of fixed effects, *a* is the vector of random breeding values, *m* is the vector of random microbial values, and *e* is the vector of random residuals. *X*, *Z*₁ and *Z*₂ are the incidence matrices for *b*, *a* and *m*. The distributions assumed for the random effects were *a* ~ N(0, *A*σ²), *m* ~ N(0, *M*σ²) and *e* ~ N(0, *I*σ²), where σ², σ²_m and σ²_m are the genetic, microbiome and residual variances, respectively; *I* denotes the identity matrix; *A* is the pedigree relationship matrix based on the 588 animals with microbiota data, plus 6705 ancestors (parents from selection generations G0 to G9 of the lines, plus their ancestors in the shared original population); and *M* is the microbial relationship matrix, defined as *M* = Z₃Z₃⁻¹, where Z₃ is a matrix with dimension n × k, with n being the number of animals with microbiome information and k the number of OTU [3]. Elements of the Z₃ matrix are the standardized individual abundance
of each OTU j for individual i, according to the following equation:

$$z_{3ij} = \frac{\log(P_{ij}) - \overline{\log(P_{ij})}_j}{sd(\log(P_{ij}))_j}, \qquad (5)$$

where Pij is the relative abundance of OTU j for individual i. The fixed effects fitted in each model were pen size (5 levels), herd of birth (2 levels), sex (3 levels), and contemporary groups (CG, 109 levels, with both lines belonging to each CG). Their significance (p<0.05) on the five traits was tested in preliminary linear models. The four models were used to evaluate their goodness-of-fit regarding the microbiome effect, using model comparisons between models with and without the microbial values with one degree of freedom, i.e. Model (1) with Model (2) and Model (3) with Model (4).

In addition, to assess if covariance between traits could be due to similar microbiota compositions [5], microbial correlations (r<sub>m</sub>) between the traits were quantified. Bivariate analyses were conducted with Models (2), (3), and (4). In this case, the distributions assumed for the random terms were **a** ~ N(0, G<sub>0</sub> ⊗ **A**), **m** ~ N(0, R<sub>m</sub> ⊗ **M**), and **e** ~ N(0, R<sub>e</sub> ⊗ **I**), where G<sub>0</sub> = [σ<sup>2</sup><sub>a<sub>i</sub></sub> σ<sub>a<sub>ij</sub></sub> σ<sup>2</sup><sub>a<sub>ji</sub></sub> σ<sup>2</sup><sub>a<sub>j</sub></sub> σ<sub>a<sub>ij</sub></sub> σ<sub>a<sub>ji</sub></sub> σ<sup>2</sup><sub>m<sub>i</sub></sub> σ<sub>m<sub>ij</sub></sub> σ<sup>2</sup><sub>m<sub>ji</sub></sub> σ<sup>2</sup><sub>m<sub>j</sub></sub> σ<sub>m<sub>ji</sub></sub> σ<sub>m<sub>j</sub></sub> σ<sup>2</sup><sub>e<sub>i</sub></sub> σ<sub>e<sub>ij</sub></sub> σ<sup>2</sup><sub>e<sub>ji</sub></sub> σ<sup>2</sup><sub>e<sub>j</sub></sub> σ<sub>e<sub>ij</sub></sub> σ<sub>e<sub>ji</sub></sub> σ<sup>2</sup><sub>e<sub>j</sub></sub> σ<sub>e<sub>ij</sub></sub> σ<sub>e<sub>ji</sub></sub> σ<sup>2</sup><sub>\*</sub> are 2 × 2 symmetric (co)variance matrices of genetic effects with the previously defined genetic variances and the genetic correlation r<sub>a<sub>ij</sub></sub> = σ<sub>a<sub>ij</sub></sub>/σ<sub>a<sub>i</sub></sub> σ<sub>a<sub>j</sub></sub> between each pair of traits i and j, and similarly R<sub>m</sub> = [σ<sup>2</sup><sub>m<sub>i</sub></sub> σ<sub>m<sub>ij</sub></sub> σ<sub>m<sub>ji</sub></sub> σ<sup>2</sup><sub>m<sub>j</sub></sub> σ<sub>m<sub>ji</sub></sub> σ<sub>m<sub>j</sub></sub> σ<sup>2</sup><sub>e<sub>i</sub></sub> σ<sub>e<sub>ij</sub></sub> σ<sup>2</sup><sub>e<sub>ji</sub></sub> σ<sup>2</sup><sub>e<sub>j</sub></sub> σ<sub>e<sub>ij</sub></sub> σ<sub>e<sub>ji</sub></sub> σ<sup>2</sup><sub>e<sub>j</sub></sub> σ<sub>e<sub>ij</sub></sub> σ<sub>e<sub>ji</sub></sub> σ<sup>2</sup><sub>\*</sub> are 2 × 2 symmetric (co)variances matrices of microbial and residual effects, respectively.

All models were fitted using the Bayesian framework using the GIBBSF90 software [19]. In total, 100,000 samples were generated to obtain the posterior distributions of the parameters for each model, and a burn-in period of 15,000 samples and thinning interval of 10 were considered. The convergence was verified through visual inspection of trace sample plots. Model comparisons were based on the deviance information criterion (DIC; Spiegelhalter et al. [20]). As a further evaluation of models, correlations of phenotypes adjusted for fixed effects (y<sup>*</sup>) with the sum of solutions obtained for random terms (except residuals) for each individual were calculated and presented as r(y<sup>*</sup>, m̂) for Model (2), r(y<sup>*</sup>, â) for Model (3) and r(y<sup>*</sup>, â + m̂) for Model (4).

## Accuracy of microbiome predictions of phenotypes

To evaluate the accuracy of prediction of the phenotypes using m̂, two scenarios for cross-validation were considered for each of the five traits using Models (2) and (4): the first scenario was designed to run predictions for 50 random animals in 20 successive replicates, and the second scenario was designed to run predictions for the animals of each of the 14 CG, successively. The objective with these designs was to evaluate the effect of the presence of contemporary animals on the prediction accuracy.

The part-whole linear regression (LR) method of Legarra and Reverter [21] was used to quantify prediction accuracies. This is an alternative to the conventional calculation of accuracy based on correlations with adjusted phenotypes, which can be affected by the limited accuracy of the computation of adjusted phenotypes. For the LR method, the phenotypes of the animals to be predicted were set to missing for each trait and their m̂ were predicted (m̂<sub>p</sub>, i.e. prediction using a partial dataset) with Models (2) and (4). For each replicate, prediction accuracy for each trait was evaluated based on the correlation between m̂<sub>p</sub> and predicted m̂ using the full dataset (m̂<sub>w</sub>, i.e., prediction using the whole dataset) with the same model. The final criterion was the average of the correlations across replicates (r(m̂<sub>p</sub>, m̂<sub>w</sub>)). Similarly, to evaluate the accuracy of microbiome predictions of phenotypes with Model (4), the average of part-whole correlations for the sum of predicted breeding values and microbial value was calculated, i.e. r((â + m̂)<sub>p</sub>, (â + m̂)<sub>w</sub>).

## Microbiome-wide association studies (MWAS)

The objective of the MWAS was to identify OTU that have significant associations with the phenotype of each of the five traits. Two approaches were used.

### Single-OTU regression analysis

Single-OTU regression analyses were applied to test the effect of the 2630 OTU, one at a time, and obtain the associated p-value, which is the most common approach [10]. The model used for these analyses was the same as Model (3) except that a specific OTU was added as a fixed covariate. The model was fitted using the best linear unbiased prediction (BLUP) method of the AIREMLF90 software [19]. The p-value of the estimate of the regression coefficient for the fitted OTU covariate was obtained by converting the estimate and its standard error to a Z-score and applying a Chi-squared test.
## Back-solving BLUP solutions

An alternative approach is the back-solving of predicted $\hat{\mathbf{m}}$ from Model (4) to obtain the estimated effect of each OTU on the trait phenotypes. This approach is similar to those described by Stranden and Garrick [22] and Gualdron-Duarte et al. [23] to obtain estimates and prediction error variances for single nucleotide polymorphisms (SNPs) from estimated breeding values using genomic-BLUP. Such back-solving is often used for SNP genome-wide association studies. It has only recently been proposed for microbiota analyses [11], however, the corresponding p-values of the estimates have not been reported. The back-solving approach and calculation of the p-values were conducted as described in the following.

Estimates for the effect of each OTU ($\widehat{\text{OTU}}$) can be obtained if the assumptions of $\sigma_{\text{OTU}}^2 = \sigma_m^2/k$ and $\mathbf{D} = \mathbf{I}\sigma_m^2/k$ hold, and thus:

$$Z_3DZ'_3 = M\sigma_m^2. \qquad (6)$$

Therefore, estimates for the effects of OTU given the $\hat{\mathbf{m}}$ can be achieved as follows [22]:

$$E(\widehat{\text{OTU}}) = \widehat{\text{OTU}} = DZ'_3(Z_3DZ'_3)^{-1}\hat{\mathbf{m}} = \frac{1}{k}Z_3'M^{-1}\hat{\mathbf{m}}. \qquad (7)$$

And the variance of the OTU estimates is defined as in [23]:

$$\text{var}(\widehat{\text{OTU}}) = \text{var}\left(\frac{1}{k}Z_3'M^{-1}\hat{\mathbf{m}}\right) = \frac{1}{k}Z_3'M^{-1}\text{var}(\hat{\mathbf{m}})M^{-1}Z_3\frac{1}{k}. \qquad (8)$$

The predictor error variance (PEV) of $\hat{\mathbf{m}}$ is equal to:

$$\text{PEV}(\hat{\mathbf{m}}) = \text{var}(\mathbf{m} - \hat{\mathbf{m}}) = \text{var}(\mathbf{m}) - \text{var}(\hat{\mathbf{m}}) = C^{\mathbf{mm}}\sigma_e^2 \qquad (9)$$

Therefore,

$$\text{var}(\hat{\mathbf{m}}) = \text{var}(\mathbf{m}) - C^{\mathbf{mm}}\sigma_e^2 = \mathbf{M}\sigma_m^2 - C^{\mathbf{mm}}\sigma_e^2, \qquad (10)$$

where $C^{\mathbf{mm}}$ are the diagonal elements of the sub-matrix corresponding to the random microbial values from the inverse of coefficient matrix of the mixed model equations.

Finally:

$$\text{var}(\widehat{\text{OTU}}) = \frac{1}{k}Z_3'M^{-1}(\mathbf{M}\sigma_m^2 - C^{\mathbf{mm}}\sigma_e^2)M^{-1}Z_3\frac{1}{k}. \qquad (11)$$

Then, the Z-score for the estimate for each OTU j can be obtained from the diagonal elements of $\text{var}(\widehat{\text{OTU}})$ as:

$$Z - \text{score}_j = \frac{\widehat{\text{OTU}}_j}{\sqrt{\text{var}(\widehat{\text{OTU}}_j)}}. \qquad (12)$$

The corresponding p-values can then be calculated by applying a Chi-square test to these Z-scores.

The back-solving method was run using a local script to construct and solve the mixed model equations based on the variance component estimates of Model (4) for each trait. It should be noted that both the single-OTU regression and the back-solving approach accounted for additive genetic effects in order to obtain comparable results.

## Significance thresholds for MWAS

To obtain significance thresholds for the MWAS to control the family-wise type I error rate while accounting for multiple testing using non-independent variables, a principal component analysis (PCA) was applied to the correlation matrix of OTU ($Z_3'Z_{32630 \times 2630}$) to estimate the number of independent tests, as proposed by Gao et al. [24]. The PCA showed that 428 eigenvalues captured 99.5% of the variability in the correlation matrix. Based on this, thresholds for significance and suggestive significance at 5% (−log10(0.05/428)) and 10% (−log10(0.10/428)) family-wise type I error rates were used to test the significance of the effects of OTU.

## Results

### Estimation of variance components

The results of univariate analyses of the five studied traits with the four models are in Table 2, as posterior means ± posterior 95% confidence intervals of each variance component. Comparisons of the DIC values showed a consistent improvement of the fit from Model (1) to Model (4), which had the smallest DIC for all traits. This goodness-of-fit of Model (4) was confirmed by its higher $r(y^*, \hat{\mathbf{a}} + \hat{\mathbf{m}})$, compared to the correlations between $y^*$ and the predictions from Models (2) and (3) (Table 3). The posterior means for heritability were moderate for all traits and ranged from 0.31 ± 0.13 for FCR to 0.51 ± 0.10 for BFT, with no difference between estimates from Models (3) and (4). The microbiome variances obtained with Models (2) and (4) showed a substantial contribution to the phenotypic variance of feed efficiency
<table><caption>**Table 2** Posterior means (±posterior standard deviation) of variance components, heritability, and microbiability of the production traits using the four models, and the deviance information criterion (DIC) of each model and trait</caption><thead><tr><th>Trait</th><th>Model</th><th>σ<sup>2</sup><sub>g</sub></th><th>σ<sup>2</sup><sub>m</sub></th><th>σ<sup>2</sup><sub>e</sub></th><th>σ<sup>2</sup><sub>p</sub></th><th>h<sup>2</sup></th><th>m<sup>2</sup></th><th>DIC</th></tr></thead><tbody><tr><td rowspan="4">RFI</td><td>(1)</td><td>—</td><td>—</td><td>0.020 ± 0.001</td><td>0.020 ± 0.001</td><td>—</td><td>—</td><td>— 343032380251</td></tr><tr><td>(2)</td><td>—</td><td>0.002 ± 0.002</td><td>0.017 ± 0.002</td><td>0.020 ± 0.001</td><td>—</td><td>0.12 ± 0.09</td><td>— 376715009798</td></tr><tr><td>(3)</td><td>0.006 ± 0.002</td><td>—</td><td>0.014 ± 0.002</td><td>0.020 ± 0.001</td><td>0.32 ± 0.10</td><td>—</td><td>— 481592517646</td></tr><tr><td>(4)</td><td>0.006 ± 0.002</td><td>0.002 ± 0.002</td><td>0.012 ± 0.002</td><td>0.020 ± 0.001</td><td>0.30 ± 0.10</td><td>0.11 ± 0.09</td><td>— 540693245186</td></tr><tr><td rowspan="4">FCR</td><td>(1)</td><td>—</td><td>—</td><td>0.062 ± 0.004</td><td>0.062 ± 0.004</td><td>—</td><td>—</td><td>— 64257520104</td></tr><tr><td>(2)</td><td>—</td><td>0.014 ± 0.008</td><td>0.051 ± 0.007</td><td>0.065 ± 0.005</td><td>—</td><td>0.22 ± 0.11</td><td>— 78170373902</td></tr><tr><td>(3)</td><td>0.024 ± 0.010</td><td>—</td><td>0.043 ± 0.007</td><td>0.067 ± 0.005</td><td>0.35 ± 0.13</td><td>—</td><td>— 93388582161</td></tr><tr><td>(4)</td><td>0.022 ± 0.010</td><td>0.014 ± 0.008</td><td>0.032 ± 0.009</td><td>0.070 ± 0.006</td><td>0.31 ± 0.13</td><td>0.20 ± 0.11</td><td>— 122763803675</td></tr><tr><td rowspan="4">DFI</td><td>(1)</td><td>—</td><td>—</td><td>0.051 ± 0.003</td><td>0.051 ± 0.003</td><td>—</td><td>—</td><td>— 89420640342</td></tr><tr><td>(2)</td><td>—</td><td>0.003 ± 0.003</td><td>0.050 ± 0.004</td><td>0.052 ± 0.003</td><td>—</td><td>0.06 ± 0.06</td><td>— 93351683359</td></tr><tr><td>(3)</td><td>—</td><td>—</td><td>0.030 ± 0.006</td><td>0.056 ± 0.005</td><td>0.50 ± 0.13</td><td>—</td><td>— 167531111339</td></tr><tr><td>(4)</td><td>0.030 ± 0.010</td><td>0.002 ± 0.002</td><td>0.030 ± 0.007</td><td>0.060 ± 0.005</td><td>0.48 ± 0.14</td><td>0.04 ± 0.03</td><td>— 173089799612</td></tr><tr><td rowspan="4">ADG</td><td>(1)</td><td>—</td><td>—</td><td>0.0051 ± 0.0003</td><td>0.0051 ± 0.0003</td><td>—</td><td>—</td><td>— 255536305168</td></tr><tr><td>(2)</td><td>—</td><td>0.0002 ± 0.0003</td><td>0.0050 ± 0.0004</td><td>0.0051 ± 0.0003</td><td>—</td><td>0.05 ± 0.05</td><td>— 265759148148</td></tr><tr><td>(3)</td><td>0.0024 ± 0.0009</td><td>—</td><td>0.0030 ± 0.0006</td><td>0.0054 ± 0.0005</td><td>0.45 ± 0.13</td><td>—</td><td>— 440920106295</td></tr><tr><td>(4)</td><td>0.0030 ± 0.0008</td><td>0.0001 ± 0.0002</td><td>0.0030 ± 0.0006</td><td>0.0055 ± 0.0005</td><td>0.47 ± 0.12</td><td>0.03 ± 0.03</td><td>— 471061077561</td></tr><tr><td rowspan="4">BFT</td><td>(1)</td><td>—</td><td>—</td><td>8.854 ± 0.570</td><td>8.854 ± 0.570</td><td>—</td><td>—</td><td>— 525486682</td></tr><tr><td>(2)</td><td>—</td><td>0.100 ± 0.610</td><td>8.057 ± 0.680</td><td>9.055 ± 0.604</td><td>—</td><td>0.11 ± 0.06</td><td>— 575243859</td></tr><tr><td>(3)</td><td>4.754 ± 1.301</td><td>—</td><td>4.695 ± 0.872</td><td>9.450 ± 0.750</td><td>0.50 ± 0.11</td><td>—</td><td>— 999015096</td></tr><tr><td>(4)</td><td>4.980 ± 1.280</td><td>0.228 ± 0.314</td><td>4.424 ± 0.878</td><td>9.636 ± 0.760</td><td>0.51 ± 0.10</td><td>0.02 ± 0.03</td><td>— 1071710845</td></tr></tbody></table>

σ²g: genetic variance, σ²m: microbiome variance, σ²e: residual variance, σ²p: phenotypic variance, h²: heritability, m²: microbiability
RFI residual feed intake (kg/day), FCR feed conversion ratio (kg/kg), DFI daily feed intake (kg/day), ADG average daily gain (kg/day), BFT backfat thickness (mm)

<table><caption>**Table 3** Pearson correlations of estimated breeding values (â) and estimated microbial values (m̂) with the adjusted phenotypes of the production traits (y*) for the three models</caption><thead><tr><th rowspan="2">Trait</th><th colspan="3">Model</th></tr><tr><th>(2)<br>r(y*, m̂)</th><th>(3)<br>r(y*, â)</th><th>(4)<br>r(y*, a + m̂)</th></tr></thead><tbody><tr><td>RFI</td><td>0.79</td><td>0.79</td><td>0.88</td></tr><tr><td>FCR</td><td>0.81</td><td>0.79</td><td>0.88</td></tr><tr><td>DFI</td><td>0.68</td><td>0.76</td><td>0.84</td></tr><tr><td>ADG</td><td>0.61</td><td>0.75</td><td>0.80</td></tr><tr><td>BFT</td><td>0.62</td><td>0.88</td><td>0.92</td></tr></tbody></table>

RFI residual feed intake, FCR feed conversion ratio, DFI daily feed intake, ADG average daily gain, BFT backfat thickness

traits, with posterior means for m² of, respectively, 0.12 ± 0.09 and 0.11 ± 0.09 for RFI and of 0.22 ± 0.11 and 0.20 ± 0.11 for FCR. In contrast, the phenotypic variances of DFI, BFT and ADG were less influenced by the microbiome variance. The posterior means for m² were lower than 0.06 ± 0.06 for DFI and ADG for Models (2) and (4), and for BFT it ranged from 0.11 ± 0.06 with Model (2) to 0.02 ± 0.03 with Model (4), i.e. mainly close

to zero. All traits showed lower posterior means for m² than for h².

The results of the bivariate analyses for the studied traits with Models (2), (3), and (4) are in Table 4. Posterior means for h² and m² of the traits in these analyses were similar to those obtained with the univariate analyses. Given the posterior standard deviations, estimates of Models (2) and (3) were not different from the corresponding estimates of Model (4), except for the estimate for r_m between ADG and BFT. Estimates of r_m from Model (4) ranged from -0.54 ± 0.60 for RFI and ADG to 0.96 ± 0.11 for ADG and BFT. Except for the estimate of r_m between RFI and FCR (0.70 ± 0.34), all other estimates had high standard errors given the low microbiability estimates of the traits. Estimates of genetic correlations, r_a, were very similar for Models (3) and (4).

### Accuracy of microbiome predictions of phenotypes

The r(m̂p, m̂w) for 50 random animals in 20 replicates for the five traits are shown in Fig. 1 and the r(m̂p, m̂w) for CG are given in Fig. S1 (see Additional file 1: Fig. S1). The average correlations ± SD ranged from 0.55 ± 0.14 to 0.81 ± 0.05 for the random design and from 0.39 ± 0.23 to 0.70 ± 0.15 for the CG design,
<table><caption>**Table 4** Posterior means (±posterior standard deviation) of the heritability, microbiability obtained from bivariate analyses between production traits with models (2), (3) and (4)</caption><thead><tr><th>Model</th><th>Trait 1</th><th>Trait 2</th><th>h<sup>2</sup><sub>T1</sub></th><th>h<sup>2</sup><sub>T2</sub></th><th>m<sup>2</sup><sub>T1</sub></th><th>m<sup>2</sup><sub>T2</sub></th><th>r<sub>a12</sub></th><th>r<sub>m12</sub></th></tr></thead><tbody><tr><td rowspan="12">(2)</td><td rowspan="4">RFI</td><td>FCR</td><td>—</td><td>—</td><td>0.12±0.13</td><td>0.25±0.09</td><td>—</td><td>0.75±0.33</td></tr><tr><td>DFI</td><td>—</td><td>—</td><td>0.26±0.11</td><td>0.04±0.05</td><td>—</td><td>0.87±0.33</td></tr><tr><td>ADG</td><td>—</td><td>—</td><td>0.12±0.10</td><td>0.09±0.07</td><td>—</td><td>—0.87±0.39</td></tr><tr><td>BFT</td><td>—</td><td>—</td><td>0.23±0.10</td><td>0.01±0.00</td><td>—</td><td>—1.00±NE</td></tr><tr><td rowspan="4">FCR</td><td>DFI</td><td>—</td><td>—</td><td>0.23±0.12</td><td>0.12±0.07</td><td>—</td><td>0.66±0.46</td></tr><tr><td>ADG</td><td>—</td><td>—</td><td>0.17±0.11</td><td>0.13±0.07</td><td>—</td><td>—0.94±0.13</td></tr><tr><td>BFT</td><td>—</td><td>—</td><td>0.24±0.13</td><td>0.06±0.07</td><td>—</td><td>0.76±0.40</td></tr><tr><td>DFI</td><td>ADG</td><td>—</td><td>—</td><td>0.03±0.04</td><td>0.13±0.07</td><td>—</td><td>—0.38±0.66</td></tr><tr><td rowspan="2">DFI</td><td>ADG</td><td>—</td><td>—</td><td>0.06±0.04</td><td>0.16±0.07</td><td>—</td><td>0.88±0.32</td></tr><tr><td>BFT</td><td>—</td><td>—</td><td>0.04±0.05</td><td>0.14±0.06</td><td>—</td><td>—0.28±0.66</td></tr><tr><td rowspan="12">(3)</td><td rowspan="4">RFI</td><td>FCR</td><td>0.35±0.12</td><td>0.40±0.12</td><td>—</td><td>—</td><td>0.64±0.16</td><td>—</td></tr><tr><td>DFI</td><td>0.31±0.10</td><td>0.51±0.14</td><td>—</td><td>—</td><td>0.49±0.22</td><td>—</td></tr><tr><td>ADG</td><td>0.33±0.10</td><td>0.52±0.13</td><td>—</td><td>—</td><td>—0.01±0.24</td><td>—</td></tr><tr><td>BFT</td><td>0.30±0.09</td><td>0.53±0.10</td><td>—</td><td>—</td><td>0.002±0.20</td><td>—</td></tr><tr><td rowspan="4">FCR</td><td>DFI</td><td>0.33±0.15</td><td>0.52±0.14</td><td>—</td><td>—</td><td>0.50±0.26</td><td>—</td></tr><tr><td>ADG</td><td>0.39±0.13</td><td>0.47±0.14</td><td>—</td><td>—</td><td>—0.23±0.24</td><td>—</td></tr><tr><td>BFT</td><td>0.34±0.12</td><td>0.51±0.12</td><td>—</td><td>—</td><td>0.59±0.20</td><td>—</td></tr><tr><td>DFI</td><td>ADG</td><td>0.50±0.12</td><td>0.49±0.13</td><td>—</td><td>—</td><td>0.64±0.16</td><td>—</td></tr><tr><td rowspan="2">DFI</td><td>ADG</td><td>0.52±0.14</td><td>0.55±0.12</td><td>—</td><td>—</td><td>0.65±0.14</td><td>—</td></tr><tr><td>BFT</td><td>0.54±0.15</td><td>0.53±0.11</td><td>—</td><td>—</td><td>0.30±0.19</td><td>—</td></tr><tr><td rowspan="12">(4)</td><td rowspan="4">RFI</td><td>FCR</td><td>0.35±0.11</td><td>0.38±0.13</td><td>0.16±0.10</td><td>0.23±0.10</td><td>0.66±0.16</td><td>0.70±0.34</td></tr><tr><td>DFI</td><td>0.29±0.09</td><td>0.54±0.13</td><td>0.18±0.12</td><td>0.06±0.04</td><td>0.63±0.17</td><td>0.71±0.47</td></tr><tr><td>ADG</td><td>0.33±0.10</td><td>0.51±0.13</td><td>0.09±0.08</td><td>0.10±0.05</td><td>0.00±0.27</td><td>—0.54±0.60</td></tr><tr><td>BFT</td><td>0.29±0.10</td><td>0.52±0.11</td><td>0.17±0.08</td><td>0.05±0.03</td><td>0.02±0.27</td><td>—1.00±NE</td></tr><tr><td rowspan="4">FCR</td><td>DFI</td><td>0.32±0.12</td><td>0.52±0.12</td><td>0.28±0.10</td><td>0.08±0.03</td><td>0.44±0.25</td><td>0.99±NE</td></tr><tr><td>ADG</td><td>0.38±0.13</td><td>0.51±0.13</td><td>0.22±0.10</td><td>0.11±0.05</td><td>—0.25±0.21</td><td>—0.91±0.18</td></tr><tr><td>BFT</td><td>0.34±0.11</td><td>0.49±0.10</td><td>0.23±0.09</td><td>0.05±0.04</td><td>0.52±0.20</td><td>0.40±0.64</td></tr><tr><td>DFI</td><td>ADG</td><td>0.48±0.13</td><td>0.49±0.12</td><td>0.04±0.04</td><td>0.12±0.06</td><td>0.62±0.16</td><td>—0.37±0.56</td></tr><tr><td rowspan="2">DFI</td><td>ADG</td><td>0.50±0.13</td><td>0.51±0.10</td><td>0.04±0.04</td><td>0.06±0.04</td><td>0.61±0.15</td><td>0.50±0.58</td></tr><tr><td>BFT</td><td>0.49±0.14</td><td>0.50±0.11</td><td>0.04±0.05</td><td>0.06±0.05</td><td>0.30±0.20</td><td>0.96±0.11</td></tr></tbody></table>

$h_{T1}^2$: heritability of first trait, $h_{T2}^2$: heritability of second trait, $m_{T1}^2$: microbiability of first trait, $m_{T2}^2$: microbiability of second trait, $r_{a12}$: genetic correlation, $r_{m12}$: microbial correlation, NE: not estimable

RFI residual feed intake, FCR feed conversion ratio, DFI daily feed intake, ADG average daily gain, BFT backfat thickness

depending on the trait. The $\overline{r}(\hat{\mathbf{m}}_p, \hat{\mathbf{m}}_w)$ were systematically slightly higher with the random design than with the CG design for both Models (2) and (4), but differences were within the range of the SE. Finally, including the genetic effect in Model (4) did not improve prediction accuracy of $\hat{\mathbf{m}}$. With the random design, the $\overline{r}(\hat{\mathbf{m}}_p, \hat{\mathbf{m}}_w)$ for BFT was higher than for the other four traits, which had corresponding accuracies around 0.60. With the CG design, the $\overline{r}(\hat{\mathbf{m}}_p, \hat{\mathbf{m}}_w)$ was slightly lower for FCR and RFI than for ADG, BFT and DFI.

Estimates of $\overline{r}(\overline{(\hat{\mathbf{a}} + \hat{\mathbf{m}})_p}, \overline{(\hat{\mathbf{a}} + \hat{\mathbf{m}})_w})$ obtained with Model (4) ranged from $0.68 \pm 0.12$ to $0.78 \pm 0.05$ with the random design (Fig. 2) and from $0.60 \pm 0.19$ to $0.75 \pm 0.10$ with the CG design (see Additional file 1:

Fig. S2). With both designs, the highest value was obtained for ADG and the lowest value for FCR. Although not significantly different, these values were higher than the corresponding $\overline{r}(\hat{\mathbf{m}}_p, \hat{\mathbf{m}}_w)$ values for all traits, except for BFT with the random design, for which adding $\hat{\mathbf{a}}$ resulted in slightly lower accuracies (from 0.74 to 0.68).

Microbiome-wide association studies

The results of MWAS with single-OTU regression are shown in Fig. 3 (RFI and FCR) and Fig. 4 (DFI, ADG, and BFT), and those from the back-solving approach are given in Figure S3 (see Additional file 1: Fig. S3). The two approaches used for MWAS resulted in similar estimates,
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Average correlations between $\hat{m}_p$ ($\hat{m}$ using partial datasets) and $\hat{m}_w$ ($\hat{m}$ using the whole dataset) for 50 random animals in 20 replicates, and their SD as error bars. ADG average daily gain, BFT backfat thickness, DFI daily feed intake, FCR feed conversion ratio, RFI residual feed intake, EMV estimated microbiota values</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Average correlations between $(\hat{a} + \hat{m})_p$ and $(\hat{a} + \hat{m})_w$ for 50 random animals in 20 replicates, and their SD as error bars. ADG average daily gain, BFT backfat thickness, DFI daily feed intake, FCR feed conversion ratio, RFI residual feed intake</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Results of microbiome wide association analyses using single-OTU regression method between operational taxonomic units and residual feed intake (RFI) and feed conversion ratio (FCR). In the plots, the solid and dashed lines represent significance and suggestive significance at 5 and 10% family-wise type I error rates, respectively</figcaption></figure>

but the single-regression method showed slightly higher significance levels, although Rubio et al. [25] have demonstrated the theoretical equivalence between these two approaches. There were no common significant or suggestive OTU between the traits. For RFI, the single-OTU regression resulted in three suggestive OTU (OTU391, OTU1749, and OTU2280), while the back-solving method resulted in one significant (OTU391) and one suggestive OTU (OTU1749). Both approaches identified one significant OTU for FCR (OTU1768) and one for BFT (OTU2934). For DFI, the single-OTU regression identified two significant (OTU694 and OTU1619) and one suggestive OTU (OTU2678), while the back-solving method resulted in one significant (OTU694) and one suggestive OTU (OTU1619). For ADG, none of the methods detected significant or suggestive OTU.

The eight OTU with significant or suggestive associations with at least one of the five traits belonged to the *Streptococcaceae* (1 OTU), *Prevotellaceae* (3 OTU), *Ruminococcaceae* (3 OTU), and *Lachnospiraceae* (1

OTU) families (Table 5). From these, the six OTU with an identified genus belonged to different genera. All these genera had more than 85% of zero count, and only OTU391, associated with RFI, had on average more than 10 counts per sample when it was present. The estimated regression coefficients (Table 3) indicated that higher abundances of this OTU and of OTU1749 were associated with greater efficiency (reduced RFI). Higher abundances of OTU2280, OTU1768, and OTU2934 were associated with greater RFI, FCR, and BFT, respectively, while an increase in abundance in each of the detected OTU for DFI was associated with reduced feed intake.

## Discussion

### Estimates of variance components

Previous studies in pigs revealed that the abundance of some components of the microbial community are heritable, and that heritability estimates vary at different stages of pig growth [3, 10, 26], which would provide a lever to select these GIT microbial components across
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Results of microbiome wide association analyses using single-OTU regression method between operational taxonomic units and daily feed intake (DFI), average daily gain (ADG) and back fat thickness (BFT). In the plots, the solid and dashed lines represent significance and suggestive significance at 5 and 10% family-wise type I error rates, respectively</figcaption></figure>

generations. Such a stability of the microbial community across generations could be used to orientate the phenotypes of the host animals. Therefore, in the present study, gut microbial information of two divergent pig lines was fitted into linear animal mixed models to quantify the proportion of the phenotypic variance of feed efficiency

and other performance traits captured by this new component. The analyses showed that the microbial variance was substantial for traits related to feed efficiency. The estimate of m² obtained for RFI in our study was lower than the value (0.45±0.15) reported by Weishaar
<table><caption>**Table 5** Taxonomy, average counts per sample and estimates of regression coefficients from single-OTU regression for the OTU with significant or suggestive associations with the production traits</caption><thead><tr><th>OTU</th><th>Phylum</th><th>Class</th><th>Order</th><th>Family</th><th>Genus</th><th>%Zero</th><th>Average counts</th><th>Average counts for non-zero samples</th><th>Associated trait</th><th>Regression coefficient</th></tr></thead><tbody><tr><td>OTU391</td><td>Firmicutes</td><td>Bacilli</td><td>Lactobacillales</td><td>Streptococcaceae</td><td>Streptococcus</td><td>92.7</td><td>2.33</td><td>31.88</td><td>RFI</td><td>-0.09</td></tr><tr><td>OTU1749</td><td>Bacteroidetes</td><td>Bacteroidia</td><td>Bacteroidales</td><td>Prevotellaceae</td><td>Prevotella_9</td><td>98.8</td><td>0.11</td><td>9.00</td><td></td><td>-0.10</td></tr><tr><td>OTU2280</td><td>Firmicutes</td><td>Clostridia</td><td>Clostridiales</td><td>Ruminococcaceae</td><td>unknown</td><td>97.8</td><td>0.06</td><td>2.70</td><td></td><td>0.11</td></tr><tr><td>OTU1768</td><td>Firmicutes</td><td>Clostridia</td><td>Clostridiales</td><td>Ruminococcaceae</td><td>Ruminococcus I</td><td>98.5</td><td>0.14</td><td>9.22</td><td>FCR</td><td>0.17</td></tr><tr><td>OTU694</td><td>Firmicutes</td><td>Clostridia</td><td>Clostridiales</td><td>Ruminococcaceae</td><td>Ruminococcaceae UCG-014</td><td>85.0</td><td>1.00</td><td>6.67</td><td>DFI</td><td>-0.06</td></tr><tr><td>OTU1619</td><td>Bacteroidetes</td><td>Bacteroidia</td><td>Bacteroidales</td><td>Prevotellaceae</td><td>Alloprevotella</td><td>94.9</td><td>0.19</td><td>3.70</td><td></td><td>-0.14</td></tr><tr><td>OTU2678</td><td>Firmicutes</td><td>Clostridia</td><td>Clostridiales</td><td>Lachnospiraceae</td><td>XBB1006</td><td>98.3</td><td>0.03</td><td>1.90</td><td></td><td>-0.25</td></tr><tr><td>OTU2934</td><td>Bacteroidetes</td><td>Bacteroidia</td><td>Bacteroidales</td><td>Prevotellaceae</td><td>unknown</td><td>98.6</td><td>0.02</td><td>1.50</td><td>BFT</td><td>4.46</td></tr></tbody></table>

RFI residual feed intake (kg/day). FCR feed conversion ratio (kg/kg). DFI daily feed intake (kg/day). ADG average daily gain (kg/day). BFT backfat thickness (mm)
et al. [26]. For FCR, the estimates of m² were within the range of those reported by Camarinha-Silva et al. [3] (0.21 ± 0.14) and by Weishaar et al. [27] (0.13 ± 0.10). For the other performance traits, low estimates of m² were obtained, in spite of the lower DIC obtained for Models (2) and (4) than for Models (1) and (3), respectively. In their study, Camarinha-Silva et al. [3] reported a non-significant m² estimate of 0.16 ± 0.10 for feed intake, which was higher than our estimates of m² with Models (2) and (4) for DFI, but not significantly different. Camarinha-Silva et al. [3] and Weishaar et al. [27] reported moderate estimates of m² of 0.28 ± 0.13 and 0.24 ± 0.11, respectively, for ADG, which were higher than our estimates for this trait. Khanal et al. [5], in a study on the microbiability of meat quality and carcass composition traits in swine, found that the estimate of m² for backfat depth increased with age at sampling, i.e. 0.01 ± 0.02 at weaning, 0.12 ± 0.04 at mid-test, and 0.25 ± 0.04 at off-test. Our estimates of m² for BFT with the full Model (4) are comparable with their value at weaning, although the sampling time in our study is equivalent to their mid-test sampling. Compared to our study, these three previous studies used different genetic types (Piétrain or commercial crossbreds), sample sizes, and times of collection, which may explain part of the observed differences. In addition, differences in the bioinformatics processing and clustering rules of the sequences to obtain OTU tables are a source of heterogeneity between studies. The aim of our choice to work with ASV, i.e. without a clustering step, was to reduce part of this heterogeneity and to facilitate comparisons between studies [28].

To our knowledge, except for one study on meat quality and carcass composition traits [5] that considered carcass ADG and fat depth, microbial correlations between the traits studied here have not been reported in pigs. The positive high estimate of $r_m$ between RFI and FCR suggests that a common microbial community influences both traits. Khanal et al. [5] observed a decrease in some estimates of genomic correlations between meat quality and carcass traits in pigs when a microbial correlation was included in the models, and argued that the genomic correlations between traits usually estimated could be partially due to correlations between the gut microbiota composition of animals. However, given that we had previously observed significant genetic correlations between the microbial components and the studied traits [12], the reverse hypothesis may be also relevant, i.e. that part of the $r_m$ between the traits may be due to the high genetic correlations between the traits. Nevertheless, none of these phenomena were actually observed in the current study, since the posterior means of the correlations were quite stable across models, which suggests that estimates

of $r_a$ and $r_m$ between the traits do not depend on each other.

## Accuracy of microbiome predictions of phenotypes

Our objective was to investigate if inclusion of microbial information collected at mid-test in the BLUP models can improve phenotype predictions of animals based on $\hat{m}$ only (Model 2) or $\hat{a} + \hat{m}$ (Model 4). Using only the microbiome information, prediction accuracies were slightly lower than when both the microbiome and genetic information were used (higher $r(\hat{a} + \hat{m})_p, (\hat{a} + \hat{m})_w$ than $r(\hat{m}_p, \hat{m}_w)$). In details, the prediction accuracies of the microbiome part were very similar with Models 2 and 4, suggesting that adding the genetic effect in Model 4 provided additional accuracy of prediction of phenotypes by capturing new information rather than by better identifying the genetic versus microbiome information in the models. Overall, the relatively high values of the $r(\hat{a} + \hat{m})_p, (\hat{a} + \hat{m})_w$ suggest that phenotypes of animals for the studied traits can be predicted with a reasonable accuracy with this information. These results are in line with the report of Khanal et al. [8], who observed a higher predictive ability of models that included both genetic and microbiome effects for most of the traits.

In addition, lower accuracy of microbiome predictions of phenotypes were obtained with the CG designs than with the random designs. This could be due to the high impact of the contemporary environment on the microbiome composition of animals. Vigors et al. [29] reported that microbiome composition can vary between farms because of differences in management, season of the year, and sanitary status, and of the influence of the dam's diet. A preliminary non-metric multidimensional scaling analysis revealed that CG was the main driving factor of the microbiota composition in our dataset. Given that, in our study, the samples were collected on an experimental farm under a standardized management protocol, microbiota differences between CG should be limited compared to differences between farms. Therefore, including data of animals sharing the same breeding environment as the predicted animals in the training sets may provide more accurate predictions of phenotypes when using the microbiome information. In conclusion, homogenous breeding conditions between training and predicted populations, at least at the farm level, are suggested to obtain reliable microbiome predictions. This concept can also justify combining both the microbiota and genetic components in prediction models, since similar values of $r(\hat{a} + \hat{m})_p, (\hat{a} + \hat{m})_w$ were obtained between CG and
random designs: the genetic component could then stabilize the predictions across varying environments.

## Microbiome-wide association studies

As expected from previous studies [23, 25], the results of the two approaches used here to detect the association between OTU and the phenotypic traits were highly consistent but detected slightly different numbers of associated OTU. Both methods detected four OTU associated with the traits, and the single-OTU regression pointed out four suggestive OTU whereas the back-solving method identified two. This slight difference could be due to the properties of the BLUP method, which tends to shrink the effect solutions towards the mean of the population, shrinkage that can be passed on to the estimates of the OTU effects after the back-solving. However, single-effect regression and BLUP-based methods have been shown to have an equivalent power for the association studies [30], and the same conclusion seems plausible for MWAS. In addition, the detection power in GWAS is in great part due to the extent of the linkage disequilibrium between SNPs. The compositional nature of the OTU data, i.e. the fact that OTU are constrained by an arbitrary total number of sequences [31], can be assumed to have similar effects on MWAS, as it creates a correlation structure between OTU abundancies. This correlation structure would affect both approaches, and further studies are needed to understand and explain how it would impact MWAS results and to evaluate if data that are pre-transformed to break the compositional nature of the OTU abundancies would result in different outcomes.

Even with low values of microbiome variances estimated from Model (4) for DFI, ADG and BFT, the MWAS results were equivalent between the single-OTU regression and the back-solving approaches, which can be considered as an indirect confirmation of the estimated microbiome variances. Indeed, if variances deviated strongly from the actual values, biased MWAS results could be expected from the back-solving approach.

The tests indicated that the abundance of some of the OTU may be associated with the variability of phenotypic traits. We have previously shown that some microbiota genera differ between lines [12]. The aim of including the genetic covariance matrix in the model to run MWAS was to control the risk that the association of abundance of the significant OTU with phenotypes would be due to the divergent selection, since we could not conduct separate analyses for each pig line because their power would have been too limited.

In a previous study at the genera level, using the same microbiome dataset but filtered for the most abundant genera, we showed significant genetic correlations of abundance of genera from the *Lachnospiraceae*, *Ruminococcaceae*, *Prevotellaceae* and *Streptococcaceae* families with RFI, DFI, and BFT [12]. Our findings that abundances of OTU from these families are significantly associated with phenotypic traits are consistent with these previous results. However, the negative estimate of the regression coefficient for the OTU pertaining to the *Streptococcus* genus in our present MWAS is opposite to the positive genetic correlation that was estimated between RFI and the *Streptococcus* genus in our earlier study. This suggests potential antagonistic relationships between microbiota and traits at the phenotypic and genetic levels. Weishaar et al. [27] also reported that abundance of OTU from the *Lachnospiraceae* and *Prevotellaceae* families had strong effects on FCR and RFI, but without indication of the direction of these effects. The *Prevotellaceae*, *Lachnospiraceae* and *Ruminococcaceae* families are involved in the digestion of fibrous material and provide short-chain fatty acids to the host [32, 33]. Bacteria from the *Streptococcaceae* family are known to be lactic acid producer bacteria [34] that have an important role in the production of dietary enzymes, such as amylase, lipase, phytase, and protease [35]. Therefore, the OTU that were identified in the MWAS could have meaningful biological links with feed efficiency and other performance traits. We did not find any OTU that was strongly associated with more than one trait and that could be claimed as a "major OTU" affecting production trait. In addition, most of these OTU were rare, with only one, OTU391, appearing to be sufficiently abundant in our conditions to be quantified systematically with a reasonable sequencing depth. If this link between OTU391 and RFI was confirmed in more diverse conditions and at the genetic level, this OTU could be used as a biomarker in selection programs to improve feed efficiency in pigs. As a next step, it will be necessary to evaluate how accounting for the microbiota composition in linear mixed models will improve the prediction accuracies of breeding values.

## Conclusions

We have shown that microbiota information can be used to better predict traits in pigs, especially feed efficiency traits. The sizable m² and the identification of some OTU with abundances that are associated with the phenotype traits indicate that some of the microbiota components are associated with the variability of production traits. In addition, the lower accuracy of microbiome predictions of phenotypes when none of the individuals from
the contemporary group of the pigs were included in the training set suggests that connecting breeding conditions between the training and predicted datasets is needed. Altogether, prediction accuracies of phenotypes accounting for microbiome and genetic covariance between animals suggest that phenotypes of animals can be reliably predicted at mid-test. These outcomes need to be confirmed in more diverse datasets with different environmental factors that may strongly modify the microbiota composition, in order to identify the limits of the use of microbiota information for the prediction of phenotypes in pigs.

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12711-022-00717-7.

**Additional file 1: Figure S1.** Average correlations between $\hat{m}_p$ ($\hat{m}$ using a partial dataset) and $\hat{m}_w$ ($\hat{m}$ using the whole dataset) for the CG design, and their SD as error bars. **Figure S2.** Average correlations between $(\hat{a} + \hat{m})_p$ and $(\hat{a} + \hat{m})_w$ for CG designs, and their SD as error bars. **Figure S3.** Results of microbiome wide association analyses using back solving of BLUP solutions between operational taxonomic units and residual feed intake (RFI), feed conversion ratio (FCR), daily feed intake (DFI), average daily gain (ADG) and back fat thickness (BFT).

## Acknowledgements

The authors would like to thank also the experimental farm staff for data collection and breeding of the animals.

## Author contributions

AA performed the statistical analyses and drafted the paper. OZ, LC and CB performed lab work and obtained raw microbial information. YB organized data collection and husbandry of animals on the farm. AA and HG designed the analyses. HG secured funding, designed data collection and supervised the study. All authors read and approved the final manuscript.

## Funding

The study received the financial support of the French National Research Agency via the MicroFeed project, under Grant ANR-16-CE20-0003.

## Availability of data and materials

The sequences are available at the Bioproject database https://www.ncbi.nlm.nih.gov/bioproject/) with accession number PRJNA701065. Other data are available upon reasonable request to the authors.

## Declarations

### Ethics approval and consent to participate

The study was conducted in accordance with the French legislation on animal experimentation and ethics at INRAE GENESI (https://doi.org/10.15454/1.5572415481185847E12, Surgères, France), on the Rouillé farm approved by the Ministry of Higher Education, Research and Innovation under reference number 86-213-01.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

## Author details

¹GenPhySE, Université de Toulouse, INRAE, ENVT, 31320 Castanet-Tolosan,

France. ²GenESI, INRAE, 17700 Surgères, France.
Received: 22 July 2021 Accepted: 28 March 2022

Published online: 25 April 2022

## References

1. Bergamaschi M, Maltecca C, Schillebeeckx C, McNulty NP, Schwab C, Shull C, et al. Heritability and genome-wide association of swine gut microbiome features with growth and fatness parameters. Sci Rep. 2020;10:10134.

2. Verschuren LM, Calus MP, Jansman AJ, Bergsma R, Knol EF, Gilbert H, et al. Fecal microbial composition associated with variation in feed efficiency in pigs depends on diet and sex. J Anim Sci. 2018;96:1405–18.

3. Camarinha-Silva A, Maushammer M, Wellmann R, Vital M, Preuss S, Bennewitz J. Host genome influence on gut microbial composition and microbial prediction of complex traits in pigs. Genetics. 2017;206:1637–44.

4. Difford G, Lassen J, Lovendahl P, Genes and microbes, the next step in dairy cattle breeding. In: Proceedings of the 67th Annual Meeting of the European Association for Animal Production: 29 August-2 September 2016; Belfast, 2016.

5. Khanal P, Maltecca C, Schwab C, Fix J, Tiezzi F. Microbiability of meat quality and carcass composition traits in swine. J Anim Breed Genet. 2021;138:223–36.

6. Tang S, Xin Y, Ma Y, Xu X, Zhao S, Cao J. Screening of microbes associated with swine growth and fat deposition traits across the intestinal tract. Front Microbiol. 2020;11:2475.

7. Saborio-Montero A, Gutiérrez-Rivas M, López-García A, García-Rodríguez A, Atxaerandio R, Goiri I, et al. Holobiont effect accounts for more methane emission variance than the additive and microbiome effects on dairy cattle. Livest Sci. 2021;250:104538.

8. Khanal P, Maltecca C, Schwab C, Fix J, Bergamaschi M, Tiezzi F. Modeling host-microbiome interactions for the prediction of meat quality and carcass composition traits in swine. Genet Sel Evol. 2020;52:41.

9. Pérez-Enciso M, Zingaretti LM, Ramayo-Caldas Y, de Los CG. Opportunities and limits of combining microbiome and genome data for complex trait prediction. Genet Sel Evol. 2021;53:65.

10. Difford GF, Plichta DR, Løvendahl P, Lassen J, Noel SJ, Højberg O, et al. Host genetics and the rumen microbiome jointly associate with methane emissions in dairy cows. PLoS Genet. 2018;14:e1007580.

11. Vollmar S, Wellmann R, Borda-Molina D, Rodehutscord M, Camarinha-Silva A, Bennewitz J. The gut microbial architecture of efficiency traits in the domestic poultry model species Japanese quail (*Coturnix japonica*) assessed by mixed linear models. G3 (Bethesda). 2020;10:2553–62.

12. Aliakbari A, Zemb O, Billon Y, Barilly C, Ahn I, Riquet J, et al. Genetic relationships between feed efficiency and gut microbiome in pig lines selected for residual feed intake. J Anim Breed Genet. 2021;138:491–507.

13. Aliakbari A, Delpuech E, Labrune Y, Riquet J, Gilbert H. The impact of training on data from genetically-related lines on the accuracy of genomic predictions for feed efficiency traits in pigs. Genet Sel Evol. 2020;52:57.

14. Gilbert H, Bidanel J-P, Gruand J, Caritez J-C, Billon Y, Guillouet P, et al. Genetic parameters for residual feed intake in growing pigs, with emphasis on genetic relationships with carcass and meat quality traits. J Anim Sci. 2007;85:3182–8.

15. Magoč T, Salzberg SL. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics. 2011;27:2957–63.

16. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.

17. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2012;41:D590–6.

18. Rothschild D, Weissbrod O, Barkan E, Kurilshikov A, Korem T, Zeevi D, et al. Environment dominates over host genetics in shaping human gut microbiota. Nature. 2018;555:210–5.

19. Misztal I, Tsuruta S, Lourenco D, Masuda Y, Aguilar I, Legarra A. Manual for BLUPF90 family of programs. Athens: University of Georgia; 2018.
20. Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A. Bayesian measures of model complexity and fit. J R Stat Soc Series B Stat Methodol. 2002;64:583–639.

21. Legarra A, Reverter A. Semi-parametric estimates of population accuracy and bias of predictions of breeding values and future phenotypes using the LR method. *Genet Sel Evol.* 2018;50:53.

22. Strandén I, Garrick DJ. Derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit. *J Dairy Sci.* 2009;92:2971–5.

23. Gualdrón Duarte JL, Cantet RJ, Bates RO, Ernst CW, Raney NE, Steibel JP. Rapid screening for phenotype-genotype associations by linear transformations of genomic evaluations. BMC Bioinformatics. 2014;15:246.

24. Gao X, Starmer J, Martin ER. A multiple testing correction method for genetic association studies using correlated single nucleotide polymorphisms. *Genet Epidemiol.* 2008;32:361–9.

25. Bernal Rubio YL, Gualdrón Duarte JL, Bates R, Ernst C, Nonneman D, Rohrer G, et al. Meta-analysis of genome-wide association from genomic prediction models. *Anim Genet.* 2016;47:36–48.

26. Mach N, Berri M, Estellé J, Levenez F, Lemonnier G, Denis C, et al. Early-life establishment of the swine gut microbiome and impact on host phenotypes. *Environ Microbiol Rep.* 2015;7:554–69.

27. Weishaar R, Wellmann R, Camarinha-Silva A, Rodehutscord M, Bennewitz J. Selecting the hologenome to breed for an improved feed efficiency in pigs—a novel selection index. *J Anim Breed Genet.* 2020;137:14–22.

28. Callahan BJ, McMurdie PJ, Holmes SP. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. *ISME J.* 2017;11:2639–43.

29. Vigors S, Doherty JVO, Sweeney T. Colonic microbiome profiles for improved feed efficiency can be identified despite major effects of farm of origin and contemporary group in pigs. *Animal.* 2020;14:2472–80.

30. Gondro C. Populations and genetic architecture. In: Primer to analysis of genomic data using R. Cham: Springer International Publishing AG; 2015. p. 105–61. https://doi.org/10.1007/978-3-319-14475-7_4

31. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome datasets are compositional: and this is not optional. *Front Microbiol.* 2017;8:2224.

32. Biddle A, Stewart L, Blanchard J, Leschine S. Untangling the genetic basis of fibrolytic specialization by *Lachnospiraceae* and *Ruminococcaceae* in diverse gut communities. *Diversity.* 2013;5:627–40.

33. Gardiner GE, Metzler-Zebeli BU, Lawlor PG. Impact of intestinal microbiota on growth and feed efficiency in pigs: a review. *Microorganisms.* 2020;8:1886.

34. du Toit M, Huch M, Cho G, Franz CM. The genus Streptococcus. In: Holzapfel WH, Wood BJB, editors. Lactic acid bacteria: biodiversity and taxonomy. Hokoben: John Wiley & Sons; 2014. p. 457–505.

35. Kim EY, Kim YH, Rhee MH, Song JC, Lee KW, Kim KS, et al. Selection of *Lactobacillus* sp. PSC101 that produces active dietary enzymes such as amylase, lipase, phytase and protease in pigs. *J Gen Appl Microbiol.* 2007;53:111–7.

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

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