Mora et al. Genetics Selection Evolution (2022) 54:81 https://doi.org/10.1186/s12711-022-00770-2 Ge n e t i c s Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access Disentangling the causal relationship between rabbit growth and cecal microbiota through structural equation models Mónica Mora1* , María Velasco‑Galilea1,2, Juan Pablo Sánchez1, Yuliaxis Ramayo‑Caldas1 and Miriam Piles1 Abstract Background: The effect of the cecal microbiome on growth of rabbits that were fed under different regimes has been studied previously. However, the term “effect” carries a causal meaning that can be confounded because of potential genetic associations between the microbiome and production traits. Structural equation models (SEM) can help disentangle such a complex interplay by decomposing the effect on a production trait into direct host genet‑ ics effects and indirect host genetic effects that are exerted through microbiota effects. These indirect effects can be estimated via structural coefficients that measure the effect of the microbiota on growth while the effects of the host genetics are kept constant. In this study, we applied the SEM approach to infer causal relationships between the cecal microbiota and growth of rabbits fed under ad libitum ­(ADGAL) or restricted feeding ­(ADGR). Results: We identified structural coefficients that are statistically different from 0 for 138 of the 946 operational taxonomic units (OTU) analyzed. However, only 15 and 38 of these 138 OTU had an effect greater than 0.2 pheno‑ typic standard deviations (SD) on ­ADGAL and ­ADGR, respectively. Many of these OTU had a negative effect on both traits. The largest effects on ­ADGR were exerted by an OTU that is taxonomically assigned to the Desulfovibrio genus (− 1.929 g/d, CSS-normalized OTU units) and by an OTU that belongs to the Ruminococcaceae family (1.859 g/d, CSS-normalized OTU units). For ­ADGAL, the largest effect was from OTU that belong to the S24-7 family (− 1.907 g/d, CSS-normalized OTU units). In general, OTU that had a substantial effect had low to moderate estimates of heritability. Conclusions: Disentangling how direct and indirect effects act on production traits is relevant to fully describe the processes of mediation but also to understand how these traits change before considering the application of an external intervention aimed at changing a given microbial composition by blocking/promoting the presence of a particular microorganism. Background under a restricted feeding regime (­ADGR) since vari- Gut microbiota influence the metabolism of the host ation in growth rate under restricted feeding is directly and, therefore, its growth and feed efficiency (FE) [1–3], related to variation in FE because feed intake (FI) is con- which are key components of profitability, productiv- stant [1, 4]. On the other hand, average daily gain under ity, and sustainability in the meat production industry. ad libitum feeding (­ ADGAL) reflects an animal’s ability to Feed efficiency can be measured as average daily gain grow. Both, ­ADGAL and ­ADGR are heritable traits, which implies that host genetics plays an important role in reg- ulating growth and feed utilization [5]. Thus, these traits *Correspondence: monica.mora@irta.cat have been successfully used as selection criteria in exper- 1 Institute of Agrifood Research and Technology (IRTA)-Animal Breeding iments carried out in pigs and rabbits [4, 6]. In addition, and Genetics, Caldes de Montbui, Barcelona, Spain Full list of author information is available at the end of the article a genome-wide association study (GWAS) on rabbits © The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/. The Creative Commons Public Domain Dedication waiver (http://​creat​iveco​ mmons.​org/​publi​cdoma​in/​zero/1.​0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data. Mora et al. Genetics Selection Evolution (2022) 54:81 Page 2 of 14 conducted by Sánchez et al. [7] identified different quan- insights in the relationships and the biological mecha- titative trait loci (QTL) for both these traits, supporting nisms among traits than MTAM. the idea that the genetic control of the growth of animals The objectives of this research were: (1) to identify the fed under a restricted regime differs from that of ani- operational taxonomic units (OTU) with abundances in mals fed ad libitum. This result is also supported by the the fecal microbiome that impact rabbit growth ­(ADGAL estimate of the genetic correlation of 0.59 between these and ­ADGR), (2) to quantify the relative importance of traits that was reported by Piles and Sánchez [5]. In line direct and microbial-mediated host genetic effects on with results of previous research in human [8] and live- ­ADGAL and A ­ DGR, (3) to estimate the heritability of the stock [9–11] populations, a recent study conducted by abundance of the most important OTU, and (4) to assess Velasco-Galilea et al. [12] in the same rabbit population the importance of direct genetic relative to total genetic as used in [7] suggested that, in rabbit, cecal microbial covariances of ­ADGAL and ­ ADGR with abundance of composition and diversity are heritable traits, and identi- OTU. fied several host genomic regions that affect cecal micro- bial composition. However, since most of these studies are based on associations or correlations, causal relation- Methods ships between the microbiome and the traits cannot be Animals and phenotypic data implied. Thus, which fractions of the genetic effect on a The 412 animals included in the present study belong to phenotypic trait (such as ­ADGAL and ­ADGR) are exerted a line of rabbits that has been selected for post-weaning directly versus indirectly, e.g. by promoting a specific growth since 1983 and that is commonly used as a ter- microbial profile in the gut that has a favourable or unfa- minal sire line in the three-way crossbreeding scheme vourable effect on the trait, remains unknown. As Gar- for rabbit meat production [14]. They were randomly diner et al. [3] pointed out in their review: “some specific selected from five batches of a larger experiment [5]. examples are outlined where an increased abundance of Most were produced in four batches in a semi-open-air certain microbiota members may be a result of improved facility of the Institute of Agrifood Research and Tech- productivity traits rather than the cause of them”. In addi- nology (IRTA, Barcelona Spain) during the first semes- tion, when analysing the impact of the microbiome on ter of 2014, and the remaining were produced in a single a host phenotype, this effect can be confounded with batch in another facility of IRTA under controlled envi- host genetics effects if the genetic association between ronmental conditions in the spring of 2016. The animals the microbiome and the phenotype is not considered in bred in the first facility were housed in groups of eight the model. How important these indirect host genetics kits from weaning (32 days of age) until the end of the effects that are exerted through the gut microbiota are fattening period (66 days of age). The kits raised in the relative to the direct genetic effects on growth and FE has second facility were housed in groups of six kits and their not been evaluated in rabbits. Direct and indirect effects growing period was slightly shorter (from 32 to 60 days could be opposite in sign, leading to a null effect of the of age). A maximum of two kits from the same litter were genotype on the trait if they are of similar magnitude. assigned to the same cage to avoid confounding between This could be the case of a microorganism that has a neg- cage and maternal effects. Apart from the afore-men- ative effect on a phenotypic trait and for which its rela- tioned differences, all animals were raised under the same tive abundance is associated with a host genetic marker management conditions and were fed with a standard that is linked to a gene that positively affects the same pelleted diet supplemented with antibiotics, except for 23 trait through some metabolic process. If both effects are kits raised in 2016 that received the same food but free of similar magnitude, the effect of that marker would not of antibiotics. Water was provided ad libitum and feed be captured by a standard GWAS for the trait. Structural was supplied once per day in a feeder with three places equation models (SEM) [13] can disentangle the com- during the 4- to 5-week fattening period. At weaning, the plex interplay that exists between the host genome and animals were randomly assigned to the ad libitum (AL) its microbiome. SEM can also help to predict the effect or restricted (R) feeding regime. The amount of feed sup- that an external intervention would have on the system, plied to the animals under R during each week for each such as one directed at inhibiting or blocking the effect batch was computed as 0.75 times the average FI of kits of certain microorganisms. Multiple-trait animal mod- under AL from the same batch during the previous week, els (MTAM) capture correlations or associations among plus 10% to account for the increase in feed intake as the traits but do not inform about their causal relationship. animals grow [15]. To generate homogeneous groups, In contrast, SEM allow the decomposition of the total kits under each feeding regime were categorized into effect (genetic and environmental) on a trait into direct two groups according to their individual weaning weight and indirect contributions. Thus, SEM can provide better (WW) being greater or smaller than 700 g. Mora et al. Genetics Selection Evolution (2022) 54:81 Page 3 of 14 Table 1 Descriptive statistics of average daily gain under allele frequency (MAF) > 0.05; (iv) and only autosomal ad libitum ­(ADGAL) and restricted ­(ADGR) feeding SNPs with known positions in the OryCun2.0 assembly Trait N Mean SD IQR [24] were used. The final dataset contained genotypes for 114,604 SNPs on 412 rabbits. ADGAL (g/day) 218 55.09 5.91 7.57 ADGR (g/day) 194 38.84 5.27 6.89 N number of animals, SD standard deviation, IQR interquartile range Statistical analyses With the main objective of estimating the effect of the Individual body weight (BW) was recorded weekly for microbiota on growth, while taking into account that all animals. Individual ADG was computed as the slope these two traits are genetically correlated, the SEM of the within-animal regression of BW on age. This trait represented in Fig. 1 was implemented to assess the was computed for individuals under each feeding regime, causal relationship of each OTU (M) with ­ADGAL and thus obtaining ADG under AL (­ADGAL) and R (­ADGR). ­ADGR. In this model, the host genotype directly affects The number of animals and descriptive statistics for both the phenotypes (G → ­ADGAL and G → ­ADGR) but traits are in Table 1. also the gut microbiome (G → M), while, the micro- biota also affects the phenotypes (G → M → ­ADGAL Cecal sampling, microbial DNA extraction, and G → M → ­ADGR). In Fig. 1, AL←M and R←M are and bioinformatic processing the structural coefficients and indicate the effect of M For each animal, a cecal sample was collected in a sterile on ­ADGAL and A ­ DGR, respectively, while keeping G tube immediately after slaughter, kept cold in the labora- constant. Environmental factors (E) affect M and the tory (4 °C) and, stored at – 80 °C until DNA extraction. phenotypes. Full details of the DNA extraction, amplification, library For each OTU, the SEM was implemented using a preparation and sequencing are given in previous reports Bayesian analysis of three phenotypes: A ­ DGAL, ­ADGR, [16, 17]. Briefly, DNA integrity/purity was checked and OTU (referred as M). The analysis was repeated for according to the protocol of Desjardins and Conklin [18], each OTU, i.e. 946 times. The SEM for A ­ DGAL, ­ADGR and then DNA was amplified following Parada et al. [19]. ­ th OTU (m = 1, 2, …, 946) for the n animals can and the m Raw reads were processed with the QIIME software be written as: version 1.9.0 (https://​github.​com/​bioco​re/​qiime/​relea​ y = ( ⊗ In )y + Xb + Zl l + Zc c + Zu u + e, (1) ses/​tag/1.​9.0) [20], as described in [16]. Contigs with a quality score smaller than Q19 were removed. The where y is the vector of phenotypic records contain- UCHIME algorithm [21] was used to detect and remove ing observations for ­ADGAL, ­ADGR, and M (with ­ADGR chimeric sequences. Filtered sequences were clustered missing for animals on AL and ­ ADGAL missing for into OTU and those that were detected in less than 5% of the samples and with a count sum less than 0.01% were discarded. The remaining 946 OTU were normal- ized using the cumulative sum scaling (CSS) method [22]. Taxonomic assignment of each OTU was obtained by mapping the sequences against the Greengenes refer- ence database. Raw sequence data were deposited in the sequence read archive of NCBI under BioProject acces- sion number PRJNA524130. Host DNA extraction and single nucleotide polymorhism (SNP) genotyping Using the MN Nucleospin Tissue kit (Macherey-Nagel, Germany), host DNA was extracted from liver samples collected at slaughter. DNA integrity and purity were measured following the protocol of Desjardins and Conk- Fig. 1 Graphical representation of the structural equations model. lin [18]. The Affymetrix Axiom OrcunSNP array (199,692 ­ADGAL: average daily gain of growing rabbits fed ad libitum, ­ADGR: SNPs) was used to genotype 412 rabbits. Quality control average daily gain of growing rabbits under restricted feeding, M: of the SNPs was performed with the PLINK software microbiome, G: host genetics, E: environmental factors, AL←M: (version 1.9) [23] according to four criteria: (i) individual structural coefficient of the effect of M on ­ADGAL and R←M: structural coefficient of the effect of M on ­ADGR call rate > 0.90; (ii) SNP call rate > 0.95; (iii) SNP minor Mora et al. Genetics Selection Evolution (2022) 54:81 Page 4 of 14 animals on R);  is the matrix of structural coefficients, y = (I3n −  ⊗ In )−1 Xb + (I3n −  ⊗ In )−1 Zl l defined as: + (I3n −  ⊗ In )−1 Zc c + (I3n −  ⊗ In )−1 Zu u   0 0 AL←M + (I3n −  ⊗ In )−1 e  =  0 0 R←M , (2) = Xb∗ + Zl l∗ + Zc c∗ + Zu u∗ + e∗ , 0 0 0 (5) where AL←M and R←M are the effect of M on ­ADGAL which corresponds to a MTAM with: and ­ADGR, respectively; b is the vector of systematic l∗ 0 L0∗ ⊗ In 0 0 0       effects, with incidence matrix X, which included the    c∗ C0∗ ⊗ In     0    0 0 0   effects of the combination of batch and farm (5 levels)  ∼ N  ,   ,      G∗0 ⊗ G  ∗ u 0  0 0 0 and of bodyweight at weaning (2 levels: large and small); l      e∗ 0 0 0 0 R0∗ ⊗ In     and c are the vectors of permanent environmental effects (6) of litter (195 levels) and cage (189 levels), respectively; u and is the vector of genotypic effects for the 412 animals, and e is the vector of model residuals. The vectors l , c, u, and ′ L0∗ = (I3 − )−1 L0 (I3 − )−1 , (7) e are assumed to follow a joint distribution: l 0 L0 ⊗In 0 0 0       ′ c    0   0 C0 ⊗ In 0 0   C0∗ = (I3 − )−1 C0 (I3 − )−1 , (8) ∼ N , ,  u 0  0 0 G0 ⊗ G 0      e 0 0 0 0 R 0 ⊗ In ′ (3) G∗0 = (I3 − )−1 G0 (I3 − )−1 , (9) where G is the genomic relationship matrix and L0, C0, and G0 , and R0 are symmetric covariance matrices for litter, cage, genomic and residual effects, respectively, with ′ R0∗ = (I3 − )−1 R0 (I3 − )−1 . (10) elements equal to (only the upper triangular matrix is shown): In the SEM, uAL , uR, and uM (component vectors of  2  u) represent the genomic effects that directly affect σl,AL σl,AL,R σl,AL,M ­ADGAL, ­ADGR, and M, respectively. They can also be 2 L0 =  σl,R σl,R,M , described as the effect of the genome on a trait, while 2 σl,M holding the value for the remaining traits constant, i.e. the direct effects of the genome on a trait free 2   from genetic effects that are mediated through other σc,AL σc,AL,R σc,AL,M 2 phenotypic traits that have a causal influence on the C0 =  σc,R σc,R,M , 2 target trait [26]. In turn, AL←M uM and R←M uM rep- σc,M resent the indirect genomic effects that are mediated by M on ­ADGAL and ­ADGR, respectively, such that AL←M uM + uAL and uR = R←M uM + uR are 2  uAL∗ = ∗ σu,AL σu,AL,R σu,AL,M G0 =  2 σu,R σu,R,M , the overall genetic effects exerted by the genome of an 2 σu,M individual on A­ DGAL and A ­ DGR, respectively, through  2 all the paths. Thus, the indirect effects on A ­ DGAL and ­ADGR correspond to the direct effects of genes exerted  σe,AL 0 0 and R0 =  2 σe,R 0 . Note that, the parameters over M, which, in turn, influence ­ADGAL and A ­ DGR. 2 σe,M Components of G∗0 can be written as: of the SEM are identifiable for any acyclic structure σu2 ∗ ,AL = σu,AL 2 2 + 2 × AL←M × σu,M,AL + AL←M 2 , × σu,M among traits, because R0 is diagonal [25]. Solving Eq. (1), the SEM becomes: (11) σu2 ∗ ,R = σu,R 2 2 + 2 × R←M × σu,M,R + R←M 2 × σu,M , (I3n −  ⊗ In )y = Xb + Zl l + Zc c + Zu u + e, (4) (12) where I3n is a (3n × 3n) identity matrix. σu2 ∗ ,M = σu,M 2 , (13) Thus, the reduced model is: Mora et al. Genetics Selection Evolution (2022) 54:81 Page 5 of 14 ­ADGR, and 14 affecting both traits. Note that the resolu- where σu2 ∗ ,AL, σu2 ∗ ,R , and σu2 ∗ ,M are the total genetic vari- tion of the 16S rDNA locus that was used in this study to ance of A ­ DGAL, ­ADGR, and M, respectively; σu,AL 2 , σu,R 2 , characterize the cecal microbiota allowed annotation of and u,M are the genetic variance of direct effects on σ 2 only about half of all sequences at the family level. Thus, ­ADGAL, ­ADGR, and M, respectively; and σu,AL,M and most of the OTU with a positive effect on ADG could σu,R,M are the covariances between direct effects, which not be assigned at the family level (57, 42, and 50% for represent the effects of genes that directly affect both ­ADGAL, ­ADGR, and both traits, respectively). In the three traits (i.e. pleitropic A ­ DGAL and ­ ADGR effects) or of cases, the majority of these OTU belong to the Lachno- genes that affect one of the two traits and are in linkage spiraceae and Ruminococcaceae families. For these two disequilibrium with each other, respectively. Structural families, 22 and 17%, respectively, of the total number coefficients AL←M and R←M were as defined above. of OTU had a positive effect on ­ADGAL, 22 and 22% on And ­ADGR, and 29 and 14% on both traits (Fig. 2). σu∗ ,AL,R = σu,AL,R + AL←M × σu,M,R + R←M Of the OTU with negative effects, 13 affected A ­ DGAL, 2 × σu,M,AL + R←M × AL←M × σu,M , 11 affected ­ADGR, and 10 affected both traits, of which 31, 24, and 30%, respectively, were not assigned at the (14) family level. Most of these OTU are encompassed by 2 σu∗ ,AL,M = σu,AL,M + AL←M × σu,M , (15) the Lachnospiraceae (18%), Erysipelotricaceae (18%) and Ruminococcaceae families (17%) for ­ADGAL, by the and Ruminococcaceae (19%) and Erysipelotricaceae families 2 (19%) for A ­ DGR, and by the Erysipelotricaceae family σu∗ ,R,M = σu,R,M + R←M × σu,M , (16) (40%) for both traits (Fig. 3). are the total covariance between A ­ DGAL and A ­ DGR, The marginal posterior means of the positive structural between ­ADGAL and M, and between ­ ADGR and M, coefficients for the effect of M (i.e. the specific OTU) on respectively. Note that if AL←M × σu,M 2 and  2 R←M × σu,M ­ADGAL ranged from 0.369 to 1.271 (g/d CSS-normalized have opposite signs but equal magnitude to σu,AL,M OTU units) (see Additional file 1: Table S1), and from and σu,R,M, respectively, the total genetic covariances 0.452 to 1.859 (g/d CSS-normalized OTU units) for the between the ADG traits and M could be null. However, effect of M on A ­ DGR (see Additional file 1: Table S2). In σu,AL,R + AL←M × σu,M,R + R←M × σu,M,AL + R←M × AL←M turn, the marginal posterior means of the negative struc- 2 ×σu,M is the indirect genetic association between A­ DGAL tural coefficients for the effect of M on ­ADGAL ranged and ­ADGR (i.e., the part of the covariance that is due to from − 0.770 to − 1.907 (g/d CSS-normalized OTU microbial-mediated effects). The same Eqs. (14), (15), and units) (see Additional file 1: Table S1), and from − 0.707 (16) for the partition of the total genetic variances and to − 1.929 (g/d CSS-normalized OTU units) for the covariances also apply to the other variance components effect of M on A­ DGR (see Additional file 1: Table S2). To of the model. identify and compare OTU with a substantial effect on The SEM was implemented under a Bayesian approach ­ADGAL and ­ ADGR, the structural coefficient estimate via Gibbs sampling using the Gibbsf90 software of the for each OTU and trait was expressed as the expected blupf90 family of programs [27]. Single sampling pro- change in number of phenotypic standard deviations cesses of 1,000,000 iterations were run for all the models, (SD) for ADG for a one-unit increase in CSS-normalized discarding the first 300,000 iterations of each chain and OTU units by dividing the structural coefficient by the saving 1 of every 100 samples. Posterior marginal infer- SD of the trait (Fig. 4). ences on variance components were performed with the Given the large number of OTU with structural coef- Postgibbsf90 software [27]. Structural coefficients were ficients that were statistically different from 0, only considered to be statistically different from 0 when the those with an absolute effect equal or higher than 0.2 high posterior density interval at 95% of confidence level SD units are described in detail. Estimates of these (HPD95%) did not include the 0. OTU, which we will refer to as relevant OTU, are listed in Table 2 (for ­ADGAL) and 3 (for A ­ DGR). Among the OTU with a negative effect on the production traits, Results the New.ReferenceOTU4438, a member of the S24-7 Recursive effects family, was the OTU with the largest effect on A ­ DGAL Our SEM approach allowed the identification of 138 (− 1.907 CSS-normalized OTU units) (Table 2) and OTU with structural coefficients that were statistically the New.ReferenceOTU369, taxonomically assigned to different from 0: 82 OTU had an effect on ­ADGAL, 80 on the Desulfovibrio genus, was the OTU with the largest ­ADGR, and 24 on both traits. Of these 138 OTU, 104 had a positive effect, with 45 affecting ­ADGAL, 45 affecting Mora et al. Genetics Selection Evolution (2022) 54:81 Page 6 of 14 Fig. 2 a Venn diagram showing the overlap between OTU with an abundance that had a positive effect on average daily gain of growing rabbits fed ad libitum (­ ADGAL) and under restricted feeding ­(ADGR) and b distribution of the OTU, grouped by families, with a positive effect on A ­ DGAL, on ­ADGR, or on both traits Fig. 3 a Venn diagram showing the overlap between the OTU with a an abundance that had a negative effect on average daily gain of growing rabbits fed ad libitum ­(ADGAL) and under restricted feeding ­(ADGR) and b distribution of the OTU, grouped by families, with a negative effect on ­ADGAL, on ­ADGR and on both traits effect on ­ADGR (− 1.929, g/d, CSS-normalized OTU an estimated increment of 1.8 g/d per CSS-normalized units) (Table 3). Among the OTU with a positive effect OTU unit. Of the relevant OTU with a positive effect on ­ADGR, a member of the Ruminococcaceaae family on ­ADGAL, OTU 641783 that had the largest effect and (New.ReferenceOTU1337) had the largest effect, with also belonged to the Ruminococcaceaae family, with an Mora et al. Genetics Selection Evolution (2022) 54:81 Page 7 of 14 Fig. 4 Histogram of the estimates of the effect of abundance of OTU on average daily gain of rabbits fed ad libitum (­ ADGAL) or under restricted feeding ­(ADGR) in SD units of the trait Table 2 Posterior means, posterior medians, and 95% highest posterior density intervals ­(HPD95%) of structural coefficients (in g/d, CSS-normalized OTU units) for the relevant OTU on average daily gain of rabbits fed ad libitum, along with their posterior means expressed in SD units and their assignment at the lowest taxonomic level OTU Mean Median HPD95% Effect Taxonomy New.ReferenceOTU4438 − 1.907 − 1.890 [− 3.488, − 0.393] − 0.322 Family S24-7 209947 − 1.867 − 1.844 [− 3.478, − 0.352] − 0.315 Order Clostridiales New.ReferenceOTU4568 − 1.537 − 1.546 [− 2.862, − 0.328] − 0.260 Family Lachnospiraceae 339336 − 1.440 − 1.439 [− 2.456, − 0.341] − 0.244 Order Clostridiales New.ReferenceOTU3820 − 1.438 − 1.435 [− 2.344, − 0.577] − 0.243 Order Clostridiales 860192 − 1.389 − 1.389 [− 2.355, − 0.452] − 0.235 Genus Coprobacillus 522353 − 1.355 − 1.351 [− 2.312, − 0.375] − 0.229 Genus Coprobacillus New.ReferenceOTU2572 − 1.327 − 1.322 [− 2.300, − 0.455] − 0.224 Genus Coprobacillus New.ReferenceOTU1080 − 1.310 − 1.312 [− 2.093, − 0.548] − 0.221 Genus Blaultia 641783 1.276 1.276 [0.367, 2.219] 0.215 Family Ruminococcae New.ReferenceOTU782 1.271 1.268 [0.256, 2.315] 0.215 Unassigned New.ReferenceOTU776 1.252 1.249 [0.148, 2.276] 0.212 Family Lachnospiraceae 554303 − 1.245 − 1.244 [− 2.055, − 0.476] − 0.211 Family Lachnospiraceae New.ReferenceOTU2945 1.210 1.203 [0.452, 2.030] 0.205 Unassigned New.ReferenceOTU3271 1.194 1.195 [0.055, 2.461] 0.202 Order Clostridiales estimated increment of 1.3 g/d per CSS-normalized and all of these had a negative effect on both traits: OTU unit. None of the relevant OTU had opposite New.ReferenceOTU2572 (genus Coprobacillus,) New. effects on the two traits. Of the 24 OTU with an effect ReferenceOTU3820 (order Clostridiales), New.Referen- on both traits, five (21%) were among the relevant OTU ceOTU4568 (family Lachnospiraceae), 522,353 (genus Coprobacillus) and 860,192 (genus Coprobacillus). Mora et al. Genetics Selection Evolution (2022) 54:81 Page 8 of 14 Table 3 Posterior means, posterior medians, and 95% highest posterior density intervals ­(HPD95%) of structural coefficients (in g/d, CSS-normalized OTU units) for the relevant OTU on average daily gain of rabbits fed on restriction, along with their posterior means expressed in SD units and their assignment at the lowest taxonomic level OTU Mean Median HPD95% Effect Taxonomy New.ReferenceOTU369 − 1.929 − 1.904 [− 3.224, − 0.699] − 0.366 Genus Desulfovibrio New.ReferenceOTU1337 1.859 1.871 [0.158, 3.403] 0.353 Family Ruminococcaceae New.ReferenceOTU381 1.838 1.858 [0.110, 3.435] 0.349 Family Lachnospiraceae New.ReferenceOTU1863 1.833 1.837 [0.857, 2.829] 0.348 Order Clostridiales 207340 1.793 1.791 [0.573, 3.018] 0.340 Family Mogibacteriaceae New.ReferenceOTU570 − 1.704 − 1.624 [− 3.194, − 0.475] − 0.323 Genus Phascolarctobacterium New.ReferenceOTU3941 1.611 1.598 [0.745, 2.548] 0.306 Family Christensenellaceae New.ReferenceOTU2872 1.546 1.536 [0.659, 2.449] 0.293 Order Clostridiales 213671 − 1.486 − 1.488 [− 2.559, − 0.290] − 0.282 Family Rikenellaceae New.ReferenceOTU3611 − 1.481 − 1.476 [− 2.675, − 0.233] − 0.281 Family Erysipelotrichaceae New.ReferenceOTU3526 1.447 1.443 [0.624, 2.253] 0.275 Family Lachnospiraceae 522353 − 1.438 − 1.439 [− 2.494, − 0.406] − 0.273 Genus Coprobacillus New.ReferenceOTU4568 − 1.428 − 1.408 [− 2.572, − 0.352] − 0.271 Family Lachnospiraceae 644244 − 1.423 − 1.421 [− 2.532, − 0.394] − 0.270 Order Clostridiales New.ReferenceOTU1522 1.400 1.404 [0.168, 2.561] 0.266 Family Ruminococcaceae New.ReferenceOTU3977 1.352 1.305 [0.498, 2.223] 0.257 Family Lachnospiraceae New.ReferenceOTU413 1.340 1.327 [0.275, 2.315] 0.254 Family Ruminococcaceae 348609 1.317 1.309 [0.493, 2.123] 0.250 Family Christensenellaceae 211066 − 1.314 − 1.310 [− 2.396, − 0.189] − 0.249 Family Ruminococcaceae 860192 − 1.301 − 1.297 [− 2.370, − 0.293] 0.247 Genus Coprobacillus New.ReferenceOTU4280 1.269 1.247 [0.044, 2.469] 0.241 Species Gnavus New.ReferenceOTU3816 1.216 1.217 [0.051, 2.323] 0.231 Family Ruminococcaceae New.ReferenceOTU3320 1.203 1.209 [0.376, 1.950] 0.228 Order Clostridiales New.ReferenceOTU2572 − 1.194 − 1.201 [− 2.156, − 0.246] − 0.227 Genus Coprobacillus New.ReferenceOTU3360 1.189 1.203 [0.269, 2.046] 0.226 Order Clostridiales New.ReferenceOTU3820 − 1.186 − 1.181 [− 2.297, − 0.016] − 0.225 Order Clostridiales New.ReferenceOTU1988 1.173 1.173 [0.254, 2.059] 0.223 Order Clostridiales New.ReferenceOTU362 1.165 1.140 [0.026, 2.341] 0.221 Family Ruminococcaceae New.ReferenceOTU1139 1.161 1.165 [0.123, 2.222] 0.220 Family Lachnospiraceae New.ReferenceOTU4624 − 1.145 − 1.147 [− 1.969, − 0.325] − 0.217 Family Ruminococcaceae New.ReferenceOTU669 − 1.125 − 1.135 [− 2.059, − 0.125] − 0.213 Family Ruminococcaceae 266198 1.108 1.113 [0.031, 2.120] 0.210 Order RF39 New.ReferenceOTU591 1.102 1.107 [0.174, 1.971] 0.209 Family Peptococcaceae 1108356 1.097 1.096 [0.195, 2.048] 0.208 Order RF39 New.ReferenceOTU2960 − 1.093 − 1.0.94 [− 1.780, − 0.433] − 0.207 Order Clostridiales New.ReferenceOTU4631 1.093 1.106 [0.123, 2.094] 0.207 Order Clostridiales New.ReferenceOTU1502 1.066 0.068 [0.322, 1.717] 0.202 Genus Blautia New.ReferenceOTU1728 1.055 1.067 [0.098, 2.060] 0.200 Genus Coprococcus Contribution of direct and microbial‑mediated genetic 2 × R←M × σu,M,R + R←M 2 2 × σu,M for ­ADGAL and effects to the genetic variance ­ADGR, respectively) for the most relevant OTU regard- Figures 5 and 6 show the marginal posterior dis- ing their effects on A ­ DGAL (Fig. 5) and ­ADGR (Fig. 6). tribution of the total genetic variances (σu2 ∗ ,AL In general, the distributions of the total and the direct and σu2 ∗ ,R for ­ ADGAL and ­ ADGR, respectively) and genetic variances were nearly the same because of the of the genetic variances of direct (σu,AL 2 and σu,R 2 small contribution of indirect effects. However, two for ­ADGAL and A ­DGR, respectively) and indirect relevant OTU of the Clostridiales order (New.Refer- effects (2 × AL←M × σu,M,AL + AL←M 2 2 × σu,M and enceOTU1337 and New.ReferenceOTU381 belonging Mora et al. Genetics Selection Evolution (2022) 54:81 Page 9 of 14 Fig. 5 Estimates of the total additive genetic variance (σu2∗ ,AL) for average daily gain of rabbits fed ad libitum ­(ADGAL), and of direct (σu,AL 2 ) and 2 2 indirect genetic effects (2 × AL←M × σu,M,AL + AL←M × σu,M) that contribute to this variance for the most relevant OTU to the Ruminococcaceae and Lachnospiraceae fami- contributions of their indirect effects to the total genetic lies, respectively) exerted a negative indirect effect on variance were smaller. ­ADGR (Fig. 6), while OTU 209947, also of the Clostridi- ales order, negatively affected ­ADGAL. OTU that belong Heritability estimates for ADG and the relevant OTU to other orders also affected A ­ DGAL (New.Referen- Average heritability estimates across all 946 analyses for ceOTU4438, Bacteroidales order) and A ­DGR (New. the two ADG traits were low to moderate, with posterior ReferenceOTU369, Desulfovibrionales order), but the means (SD) of 0.215 (0.104) and 0.157 (0.101) for A ­ DGAL and ­ADGR, respectively. Means (SD) of the marginal Fig. 6 Estimates of the total additive genetic variance (σu2∗ ,R) of the average daily gain of rabbits under restricted feeding ­(ADGR), and of direct 2 ) and indirect effects (2 ×  effects (σu,R 2 2 R←M × σu,M,R + R←M × σu,M) that contribute to this variance of the most relevant OTU Mora et al. Genetics Selection Evolution (2022) 54:81 Page 10 of 14 Fig. 7 Posterior means of heritabilities of the abundance of OTU that have an important effect on growth rate of rabbits fed ad libitum (­ ADGAL), on rabbits on restricted feeding ­(ADGR) or in both feeding regimes (common) posterior distributions for the heritabilities for the most almost zero indirect effects. The indirect effects of the relevant OTU that affected ­ADGAL, and ­ADGR, or both, OTU 209947, which belongs to the Clostridiales order, are shown in Fig. 7. In general, the estimates of these her- on ­ADGAL were substantial and negative, while the indi- itabilities were also low to moderate. Posterior mean her- rect effects of members of the Desulfovibrio genus (New. itabilities for the relevant OTU for A­ DGAL ranged from ReferenceOTU369) and the Lachnospiraceae family 0.083 for OTU 522353 to 0.222 for OTU 209947, both (New.ReferenceOTU381) on A ­ DGR were substantial and members of the Clostridiales order. Posterior mean herit- negative and positive, respectively. ability estimates for the relevant OTU for ­ADGR ranged from 0.046 for a member of the Phascolarctobacterium Discussion genus (New.ReferenceOTU570) to 0.234 for the New. In this study, SEM were implemented to assess causal ReferenceOTU369, which belongs to the Desulfovibrio relationships between the cecal microbiota and two genus. The most heritable families were Desulfovibrion- traits related to rabbit growth and FE. The effect that a aceae (0.23) and Methanobacteriaceae (0.24), while the specific OTU exerts on the studied phenotypes repre- least heritable families were Veillonellaceae (0.04) and sents the indirect effects (i.e. the effect of the OTU on Peptococcaceae (0.08). the trait while holding host genetics and environmental effects constant). It is measured by the structural coeffi- Direct genetic covariances between ADG and OTU cient, which quantifies the expected change in the phe- Marginal posterior distributions of the total additive notype by a one-unit increase in CSS-normalized OTU genetic covariance between A ­ DGAL and M (σu∗ ,AL,M), units of a specific OTU. When direct and indirect effects of the covariance due to direct effects (σu,AL,M), and of are of the same magnitude but have opposite signs, the the covariance due to indirect effects of M on A ­ DGAL total effect on the trait becomes null. In such a scenario, a (AL←M × σu,M2 ) are represented in Fig. 8. The same dis- GWAS can fail to identify the genomic regions that affect tributions are shown in Fig. 9 for A­ DGR (σu∗ ,R,M, σu,R,M the production trait through multiple paths. This hypoth- and R←M × σu,M2 , for total, direct, and indirect covari- esis was tested in a previous study conducted by Tiezzi ances, respectively). Total and direct additive genetic et al. [28], who evaluated how host genetics can affect fat covariances were very similar for the New.Referen- deposition in pigs by affecting gut microbiome composi- ceOTU4568, OTU 339336, and New.ReferenceOTU1863 tion, which then results in a change in the phenotype. In (all members of the Clostridiales order), due to the that study, it was demonstrated that the genes that do not Mora et al. Genetics Selection Evolution (2022) 54:81 Page 11 of 14 Fig. 8 Additive genetic covariance between average daily gain in rabbits fed ad libitum (­ ADGAL) and abundance of the most important OTU: total 2 ) additive covariance (σu∗ ,AL,M), direct effects (σu,AL,M), and indirect effects (AL←M × σu,M Fig. 9 Additive genetic covariance between average daily gain in rabbits under restricted feeding (­ ADGR) and abundance of the most important 2 ) OTU: total additive covariance (σu∗ ,R,M), direct effects (σu,R,M) and indirect effects (R←M × σu,M Mora et al. Genetics Selection Evolution (2022) 54:81 Page 12 of 14 affect the phenotype directly can be identified in GWAS propionate short-chain fatty acids (SCFA), which act as in which effects mediated by the microbiome were energy sources. Propionate is a gluconeogenic substrate included in the model. This highlights the importance of in the liver and intestine, while acetate contributes to the knowing the direct and indirect effects on a phenotype synthesis of lipids [37, 38] and has been shown to play when designing sustainable strategies to improve pro- a direct role in the regulation of appetite in mice [39]. ductivity and sustainability in livestock. To our knowl- Another OTU with an abundance that had a substantial edge, the importance of direct versus indirect effects on positive effect on FE was the New.ReferenceOTU381 growth rate through microbial diversity and composition (Lachnospiraceae family). Velasco-Galilea et al. [36] have not yet been quantified. The only related study was reported five OTU of the Lachnospiraceae family (includ- by Saborio et al. [29], who implemented a SEM approach ing the New.ReferenceOTU381) to be among the most to estimate the effect of the relative abundance of rumen informative OTU to predict rabbit growth and FE. Two microbes on methane production in dairy cattle. and three of these OTU were, respectively, positively Most of the OTU that had statistically significant struc- and negatively correlated with average daily residual feed tural coefficients and that positively affected ­ADGAL and intake in rabbits fed ad libitum. ­ADGR belong to the Lachnospiraceae or Ruminococ- Genetic variances for A ­ DGAL and ­ADGR in the SEM caceae families. Conversely, the OTU that had statisti- can be interpreted as the variance due to direct genetic cally significant structural coefficients and that negatively effects on the traits free from the genetic effects medi- affected ­ADGAL belong to a wider set of families: Rumi- ated by M that have a causal effect on them. Indirect nococcaceae, Lachnospiraceae, Methanobacteriaceae, and genetic effects can act by reducing the total genetic vari- Erysipelotrichaceae. Lachnospiraceae and Ruminococ- ance if they are strong enough and opposite in sign to the caceae are both members of the Clostridiales order. Sev- genetic covariance between M and the phenotype. This eral members of this order have previously been found to could make it difficult to assess the genetic determin- be associated with FE across different livestock species, ism of a trait affected by M based on MTAM. However, a including pigs [30, 31], chickens [32] and beef cattle [33]. SEM enables the different sources of genetic variation to Erysipelotrichaceae is a family that is associated with be disentangled. In addition, it enables prediction of the lipid metabolism and has been linked to inflammation effects of external interventions on a set of variables. In [34]. However, its increase in abundance has been associ- our study, an external intervention could involve promot- ated with the most but also with the least efficient pigs, ing or blocking the presence of some cecal microbe. depending on the sampling origin that each study chose One member of the Ruminococcaceae family (New. for gut microbiota characterization [35]. ReferenceOTU1337) and one member of the Lachno- As previously mentioned, given the large number of spiraceae family (New.Reference381) had relevant effects OTU with statistically significant structural coefficients on the total genetic variance for A ­ DGR, with indirect that were identified by our SEM approach, our discus- effects that reduced the total variance, although the esti- sion will be focused on the OTU that had a substantial mates of the structural coefficients for these OTU were effect (equal to or larger than 0.2 units of SD of the trait) positive. This was due to the negative genetic correla- on the phenotype, i.e. the relevant OTU. We identified 15 tion between ­ADGR and the abundance of these OTU. and 38 OTU with relevant effects on ­ADGAL and ­ADGR, It is worth mentioning that some members of the Lach- respectively. The estimated effects for the relevant OTU nospiraceae family are butyrate-producing bacteria that for ­ADGR were larger than those for the OTU affecting have a beneficial effect on animal gut health [40] and ­ADGAL. Thus, only 18% of the OTU that affected ­ADGAL have been previously identified to have a strong effect on had a structural coefficient larger than 0.2 SD units for FE traits in pigs [41]. ­ADGAL, while 47% of the OTU that affected ­ADGR met Estimates of total heritability for A ­ DGAL and ­ADGR this criterion. This could indicate that growth rate is are similar to those estimated in previous studies for mainly determined by host genetics whereas FE is deter- the same rabbit population, with a posterior mean of mined by both the host genetics and the gut microbiota. 0.21 and 0.08, respectively, by Piles and Sánchez [5] and This is in agreement with results reported by Velasco- with a posterior mean of 0.15 and 0.09, respectively, by Galilea et al. [36], who found that the predictive value of Velasco-Galilea et al. [36]. Conversely, estimates of direct cecal microbiota was higher for ­ADGR than for A ­ DGAL. heritability based on the SEM were slightly greater than It is worth highlighting the effects of some other OTU estimates of total heritability for both traits. This could on ­ADGAL or A ­ DGR. For instance, a member of the Phas- be due to the negative effects of host genetics exerted colarctobacterium genus (New.ReferenceOTU570) had a through the microbiota on these growth traits. In gen- substantial effect on ­ADGR (− 1.70 [− 3.194, − 0.475]). eral, estimated heritabilities for the OTU with the largest Some species of this genus produce acetate and structural coefficients were low to moderate, suggesting Mora et al. Genetics Selection Evolution (2022) 54:81 Page 13 of 14 a weak host genetic control of the rabbit cecal microbi- effect on ­ADGR, followed by abundance of one member ota, which is in agreement with results of Velasco-Gali- of the Ruminococcaceae family, which positively affected lea et al. [12], who assessed the host genetics effects on this trait. In contrast, only 33% of the OTU that had a rel- the rabbit cecal microbiota by means of Bayes factors. evant effect on A­ DGAL had a positive effect on this trait. The highest heritability [0.234 (0.091)] was estimated for New.ReferenceOTU369, which belongs to the Desulfo- Supplementary Information vibrio genus and is one of the OTU that had a relevant The online version contains supplementary material available at https://​doi.​ effect on A ­ DGR. This OTU was previously reported as org/​10.​1186/​s12711-​022-​00770-2. heritable by Velasco-Galilea et al. [12]. Additional file 1: Table S1. Posterior means, posterior medians, 95% In SEM, the genetic covariance reflects the association highest posterior density intervals ­(HPD95%) structural coefficients differ‑ between direct effects and can be considered to be due to ent from zero (in g/d, CSS-normalized OTU units) that had an effect on the genes that directly affect two traits simultaneously or to average daily gain of growing rabbits fed ad libitum and their assign‑ ment at the lowest taxonomic level. Table S2. Posterior means, posterior linkage disequilibrium between pairs of genes that each medians, 95% highest posterior density intervals ­(HPD95%) structural coef‑ affect one of the two traits [13]. However, there is a sec- ficients different from zero (in g/d, CSS-normalized OTU units) that had ond source of association between the host phenotype an effect on the average daily gain of growing rabbits under restricted feeding and their assignment at the lowest taxonomic level. and M because the host genetic effect on M also affects the host phenotype indirectly. This second source of Acknowledgements covariation could even have an opposite sign compared We would like to thank Oscar Perucho, Josep Ramon and Carmen Requena to the first source (i.e., the covariance between direct (staff of Unitat de Cunicultura, IRTA) for their contribution to data recording genetic effects). A member of the Desulfovibrio genus and animal care during the experiment. (New.ReferenceOTU369) impaired ­ADGR, and the total Author contributions indirect effects that this OTU was estimated to exert on JPS and MP conceived the experimental design. JPS, MP and MVG collected the phenotype acted by reducing the total effects. In this biological samples. MVG and MP processed the samples in the laboratory. MM and MP analyzed the data, interpreted data, prepared figures and tables, particular case, direct and total covariances were oppo- and wrote the manuscript. MVG conducted the bioinformatics processing of site in sign. Thus, the genetic correlation between this microbial sequences. JPS, YRC and MVG helped to interpret and discuss the particular OTU and A ­ DGR would be positive if only the results and reviewed the manuscript. All authors read and approved the final manuscript. effects of host genetics were considered. On the contrary, when the total effects are considered, the relevant nega- Funding tive effect of the Desulfovibrio genus on A ­ DGR can be This study was part of the Feed-a-Gene project that received funding from the European Union’s H2020 program under Grant agreement no. 633531, the captured by the negative genetic correlation between the Spanish project RTI2018-097610R-I00 and the project PID2021-128173OR-C21 OTU and growth. This result is in agreement with previ- (GENEF3). MMF is a recipient of a “Formación de Personal Investigador (FPI)” ous reports in pigs [30, 42] in which the abundances of associated with the research project RTI2018-097610R-I00. YRC is recipient of a Ramon y Cajal post-doctoral fellowship (RYC2019-027244-I) from the Spanish members of the Desulfovibrio genus were estimated to be Ministry of Science and Innovation. negatively correlated with FE. Availability of data and materials The raw sequence data were deposited in the sequence read archive of NCBI Conclusions under the BioProject accession number PRJNA524130. Our study highlights the importance of knowing the direct effects that host genetics exerts on a phenotype, Declarations as well as the indirect effects that host genetics exerts Ethics approval and consent to participate through the gut microbiota, not only to fully describe The data collected for this analysis was obtained following the national regula‑ the processes of mediation, but also to understand the tions for humane care and the use of animals in research. change in a phenotype that can result from a modifica- Consent for publication tion of the microbial composition through an external Not applicable. intervention (e.g., by blocking/promoting the presence Competing interests of a particular microorganism). This is the first study to The authors declare that they have no competing interests. evaluate the direct and indirect effects exerted by host genetics on growth. Our results show that 63% of the Author details 1 Institute of Agrifood Research and Technology (IRTA)-Animal Breeding OTU with abundances that had relevant effects on ­ADGR and Genetics, Caldes de Montbui, Barcelona, Spain. 2 Centre for Research had positive effects on this trait. Abundance of one mem- in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Cerdanyola del Vallès, ber of the Desulfovibrio genus exerted the largest negative Barcelona, Spain. Mora et al. Genetics Selection Evolution (2022) 54:81 Page 14 of 14 Received: 1 July 2022 Accepted: 28 November 2022 22. Paulson JN, Stine OC, Corrada Bravo HC, Pop M. Differential abundance analysis for microbial marker-gene surveys. Nat Methods. 2013;10:1200–2. 23. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second- generation PLINK: rising to the challenge of larger and richer datasets. GigaScience. 2015;4:7. References 24. Carneiro M, Rubin CJ, Di Palma F, Albert FW, Alföldi J, Martinez Barrio A, 1. Drouilhet L, Achard CS, Zemb O, Molette C, Gidenne T, Larzul C, et al. Direct et al. Rabbit genome analysis reveals a polygenic basis for phenotypic and correlated responses to selection in two lines of rabbits selected for change during domestication. Science. 2014;345:1074–9. feed efficiency under ad libitum and restricted feeding: I. Production traits 25. Wu X-L, Heringstad B, Gianola D. Bayesian structural equation models for and gut microbiota characteristics. J Anim Sci. 2016;94:38–48. inferring relationships between phenotypes: a review of methodology, 2. Maltecca C, Lu D, Schillebeeckx C, McNulty NP, Schwab C, Shull C, et al. identifiability, and applications. J Anim Breed Genet. 2010;127:3–15. Predicting growth and carcass traits in swine using microbiome data and 26. Valente BD, Rosa GJ, Gianola D, Wu XL, Weigel K. Is structural equation machine learning algorithms. Sci Rep. 2019;9:6574. modeling advantageous for the genetic improvement of multiple traits? 3. Gardiner GE, Metzler-Zebeli BU, Lawlor PG. Impact of intestinal micro‑ Genetics. 2013;194:561–72. biota on growth and feed efficiency in pigs: a review. Microorganisms. 27. Misztal I, Tsuruta S, Lourenco DAL, Masuda Y, Aquilar I, Legarra A, et al. 2020;8:1886. Manual for BLUPF90 family programs. Athens: University of Georgia; 2018. 4. Nguyen NH, McPhee CP, Wade CM. Responses in residual feed intake 28. Tiezzi F, Fix J, Schwab C, Shull C, Maltecca C. Gut microbiome mediates in lines of Large White pigs selected for growth rate on restricted feed‑ host genomic effects on phenotypes: a case study with fat deposition in ing (measured on ad libitum individual feeding). J Anim Breed Genet. pigs. Comput Struct Biotechnol J. 2020;19:530–44. 2005;122:264–70. 29. Saborío-Montero A, Gutiérrez-Rivas M, García-Rodríguez A, Atxaerandio 5. Piles M, Sánchez JP. Use of group records of feed intake to select for feed R, Goiri I, López de Maturana E, et al. Structural equation models to efficiency in rabbit. J Anim Breed Genet. 2019;136:474–83. disentangle the biological relationship between microbiota and complex 6. Drouilhet L, Gilbert H, Balmisse E, Ruesche J, Tircazes A, Larzul C, et al. traits: Methane production in dairy cattle as a case of study. J Anim Breed Genetic parameters for two selection criteria for feed efficiency in rabbits. Genet. 2020;137:36–48. J Anim Sci. 2013;91:3121–8. 30. Aliakbari A, Zemb O, Billon Y, Barilly C, Ahn I, Riquet J, et al. Genetic 7. Sánchez JP, Legarra A, Velasco-Galilea M, Piles M, Sánchez A, Rafel O, relationships between feed efficiency and gut microbiome in pig lines et al. Genome-wide association study for feed efficiency in collective selected for residual feed intake. J Anim Breed Genet. 2021;138:491–507. cage-raised rabbits under full and restricted feeding. Anim Genet. 31. Jiang H, Fang S, Yang H, Chen C. Identification of the relationship 2020;51:799–810. between the gut microbiome and feed efficiency in a commercial pig 8. Goodrich JK, Waters JL, Poole AC, Sutter JL, Koren O, Blekhman R, et al. cohort. J Anim Sci. 2021;99:skab045. Human genetics shape the gut microbiome. Cell. 2014;159:789–99. 32. Liu J, Stewart SN, Robinson K, Yang Q, Lyu W, Whitmore MA, et al. Linkage 9. Camarinha-Silva A, Maushammer M, Wellmann R, Vital M, Preuss S, Benne‑ between the intestinal microbiota and residual feed intake in broiler witz J. Host genome influence on gut microbial composition and micro‑ chickens. J Anim Sci Biotechnol. 2021;12:22. bial prediction of complex traits in pigs. Genetics. 2017;206:1637–44. 33. Myer PR, Smith TP, Wells JE, Kuehn LA, Freetly HC. Rumen microbiome 10. Difford GF, Plichta DR, Løvendahl P, Lassen J, Noel SJ, Højberg O, et al. from steers differing in feed efficiency. PLoS One. 2015;10:e0129174. Host genetics and the rumen microbiome jointly associate with methane 34. Kaakoush NO. Insights into the role of Erysipelotrichaceae in the human emissions in dairy cows. PLoS Genet. 2018;14:e1007580. host. Front Cell Infect Microbiol. 2015;5:84. 11. Sasson G, Ben-Shabat SK, Seroussi E, Doron-Faigenboim A, Shterzer N, 35. McCormack UM, Curião T, Metzler-Zebeli BU, Magowan E, Berry DP, Reyer Yaacoby S, et al. Heritable bovine rumen bacteria are phylogenetically H, et al. Porcine feed efficiency-associated intestinal microbiota and related and correlated with the cow’s capacity to harvest energy from its physiological traits: finding consistent cross-locational biomarkers for feed. MBio. 2017;8:e00703-17. residual feed intake. mSystems. 2019;4:e00324-18. 12. Velasco-Galilea M, Piles M, Ramayo-Caldas Y, Varona L, Sánchez JP. Use of 36. Velasco-Galilea M, Piles M, Ramayo-Caldas Y, Sanchez JP. The value of gut Bayes factors to evaluate the effects of host genetics, litter and cage on microbiota to predict feed efficiency and growth of rabbits under differ‑ the rabbit cecal microbiota. Genet Sel Evol. 2022;54:46. ent feeding regimes. Sci Rep. 2021;11:19495. 13. Gianola D, Sorensen D. Quantitative genetic models for describing 37. Bose S, Ramesh V, Locasale JW. Acetate metabolism in physiology, cancer, simultaneous and recursive relationships between phenotypes. Genetics. and beyond. Trends Cell Biol. 2019;29:695–703. 2004;167:1407–24. 38. Koh A, De Vadder F, Kovatcheva-Datchary P, Bäckhed F. From dietary fiber 14. Piles M, Blasco A, Pla M. The effect of selection for growth rate on carcass to host physiology: short-chain fatty acids as key bacterial metabolites. composition and meat characteristics of rabbits. Meat Sci. 2000;54:347–55. Cell. 2016;165:1332–45. 15. Gidenne T, Combes S, Fortun-Lamothe L. Feed intake limitation strategies 39. Frost G, Sleeth ML, Sahuri-Arisoylu M, Lizarbe B, Cerdan S, Brody L, et al. for the growing rabbit: effect on feeding behaviour, welfare, performance, The short-chain fatty acid acetate reduces appetite via a central homeo‑ digestive physiology and health: a review. Animal. 2012;6:1407–19. static mechanism. Nat Commun. 2014;5:3611. 16. Velasco-Galilea M, Guivernau M, Piles M, Viñas M, Rafel O, Sánchez A, et al. 40. Bedford A, Gong J. Implications of butyrate and its derivatives for gut Breeding farm, level of feeding and presence of antibiotics in the feed health and animal production. Anim Nutr. 2017;4:151–9. influence rabbit cecal microbiota. Anim Microbiome. 2020;2:40. 41. Weishaar R, Wellmann R, Camarinha-Silva A, Rodehutscord M, Bennewitz 17. Velasco-Galilea M, Piles M, Viñas M, Rafel O, González-Rodríguez O, Guiv‑ J. Selecting the hologenome to breed for an improved feed efficiency in ernau M, et al. Rabbit microbiota changes throughout the intestinal tract. pigs-A novel selection index. J Anim Breed Genet. 2020;137:14–22. Front Microbiol. 2018;9:2144. 42. Bergamaschi M, Tiezzi F, Howard J, Huang YJ, Gray KA, Schillebeeckx C, 18. Desjardins P, Conklin D. NanoDrop microvolume quantitation of nucleic et al. Gut microbiome composition differences among breeds impact acids. J Vis Exp. 2010;45:e2565. feed efficiency in swine. Microbiome. 2020;8:110. 19. Parada AE, Needham DM, Fuhrman JA. Every base matters: assess‑ ing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ Microbiol. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub‑ 2016;18:1403–14. lished maps and institutional affiliations. 20. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequenc‑ ing data. Nat Methods. 2010;7:335–6. 21. Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIIME improves sensitivity and speed of chimera detection. Bioinformatics. 2011;27:2194–200.