Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Ge n e t i c s https://doi.org/10.1186/s12711-022-00738-2 Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access Use of Bayes factors to evaluate the effects of host genetics, litter and cage on the rabbit cecal microbiota María Velasco‑Galilea1,2*, Miriam Piles1, Yuliaxis Ramayo‑Caldas1, Luis Varona3 and Juan Pablo Sánchez1 Abstract Background: The rabbit cecum hosts and interacts with a complex microbial ecosystem that contributes to the variation of traits of economic interest. Although the influence of host genetics on microbial diversity and specific microbial taxa has been studied in several species (e.g., humans, pigs, or cattle), it has not been investigated in rabbits. Using a Bayes factor approach, the aim of this study was to dissect the effects of host genetics, litter and cage on 984 microbial traits that are representative of the rabbit microbiota. Results: Analysis of 16S rDNA sequences of cecal microbiota from 425 rabbits resulted in the relative abundances of 29 genera, 951 operational taxonomic units (OTU), and four microbial alpha-diversity indices. Each of these microbial traits was adjusted with mixed linear and zero-inflated Poisson (ZIP) models, which all included additive genetic, litter and cage effects, and body weight at weaning and batch as systematic factors. The marginal posterior distributions of the model parameters were estimated using MCMC Bayesian procedures. The deviance information criterion (DIC) was used for model comparison regarding the statistical distribution of the data (normal or ZIP), and the Bayes factor was computed as a measure of the strength of evidence in favor of the host genetics, litter, and cage effects on micro‑ bial traits. According to DIC, all microbial traits were better adjusted with the linear model except for the OTU present in less than 10% of the animals, and for 25 of the 43 OTU with a frequency between 10 and 25%. On a global scale, the Bayes factor revealed substantial evidence in favor of the genetic control of the number of observed OTU and Shan‑ non indices. At the taxon-specific level, significant proportions of the OTU and relative abundances of genera were influenced by additive genetic, litter, and cage effects. Several members of the genera Bacteroides and Parabacteroides were strongly influenced by the host genetics and nursing environment, whereas the family S24-7 and the genus Ruminococcus were strongly influenced by cage effects. Conclusions: This study demonstrates that host genetics shapes the overall rabbit cecal microbial diversity and that a significant proportion of the taxa is influenced either by host genetics or environmental factors, such as litter and/or cage. Background The bacterial communities that inhabit the rabbit gas- trointestinal tracts play a key role in its metabolism, nutrition, and state of the immune system [1]. In the par- ticular case of this herbivorous mammal, the richest and *Correspondence: maria.velasco@cragenomica.es most diverse microbial community lies in its cecum [2]. 1 Institute of Agrifood Research and Technology (IRTA)-Animal Breeding Rabbit cecal microbial composition and diversity evolve and Genetics, Caldes de Montbui, Barcelona, Spain Full list of author information is available at the end of the article from a simple and unstable community at birth into a © 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. Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 2 of 15 complex and more stable one in the adult individuals [3]. coming from a Poisson distribution with parameter  (the Although this stability is reached at adulthood, previous mean number of observations) with probability 1 − p. studies have revealed that external factors, such as feed Therefore, the objective of our study was to dissect the composition [4, 5], level of feeding [6, 7], or the adminis- effects of the host genetics, litter and cage on a set of 984 tration of antibiotics [6–8] have a role in shaping the gut microbial traits that were defined so that they represent microbial composition and diversity. In addition to these the rabbit cecal microbiota at different levels of complex- environmental factors, host genetics could also play an ity (i.e., relative abundances of 29 genera, 951 normal- important role. Several studies in humans [9, 10], cattle ized OTU, and four microbial alpha-diversity indices) [11–15], pigs [16–19], and mice [20] have investigated in a meat rabbit population raised under standard com- the role of the host genetics on gut microbiota and have mercial conditions. These traits were analyzed using reported moderate heritabilities for certain microbial Bayesian linear and ZIP mixed models, and the statistical taxa and diversity indices. Thus, there is growing interest relevance of the ratios of the different variance compo- in investigating the interplay between the host genetics nents to the phenotypic variance estimates was evaluated and the gut microbiota and its impact on many complex through Bayes factors (BF). traits, such as human diseases, feed efficiency, or meth- ane emissions in cattle. Methods In rabbit breeding, production traits such as feed effi- Animals ciency and growth are key elements of economic profit In total, 425 meat rabbits from the Caldes line [26] were [21]. Studies that attempt to unravel the existence of a included in this study that was conducted at the Insti- potential link between these traits and the host genet- tute of Agrifood Research and Technology (IRTA). Of ics and microbiota are of great relevance for the rabbit these 425 rabbits, 336 were produced in four batches industry to define effective genetic selection and produc- and housed in collective cages, each containing eight tion strategies that contribute to sustainable production kits, in a semi-open-air facility during the first semes- and animal well-being. In this respect, previous studies ter of 2014, and 89 were produced in a single batch and have reported associations between the gut microbiota housed in collective cages, each containing six kits, in and growth [22] or feed efficiency in rabbits [23]. In addi- another facility under better controlled environmen- tion, a substantial percentage of the phenotypic variance tal conditions during the spring of 2016. After weaning of growth, feed intake and feed efficiency in growing rab- (32 days of age), all the animals were kept under the same bits has been attributed to variation in their cecal micro- management conditions and were fed with a standard biota [24]. However, to initiate selective breeding for the pelleted diet supplemented with antibiotics, except 23 presence of microbial taxa that are positively associated rabbits from the second facility that received a diet free with relevant traits, the genetic background of the rabbit of antibiotics. The fattening period lasted 5 and 4 weeks cecal microbiota needs to be dissected. Velasco-Galilea for the animals raised in the first and the second facility, et al. [24] has already provided some indirect evidence of respectively, and during the last fattening week, all the host genetics control of the rabbit cecal microbiota since animals received an antibiotic-free diet. Water was sup- the predictive value of the microbial information for feed plied ad libitum and food was provided once per day in a efficiency and other performance traits can be partially feeder with three places. After weaning, kits were classi- explained by the host additive genetic effects. Nonethe- fied into two groups according to their size (“big” if their less, it is necessary to explicitly assess whether an overall body weight was greater than 700 g or “small” otherwise) host genetics control of microbiota exists or, on the con- and randomly assigned to ad libitum (AL) feeding regime trary, whether only certain taxa or operational taxonomic or a regime restricted (R) to 75% of the AL feed intake. units (OTU) are influenced by genetic effects. Moreover, The amount of feed supplied to the animals under R for to design effective breeding programs based on microbial each week and each batch was computed as 0.75 times information, it is necessary to know whether heritable the average feed intake of kits on AL from the same batch taxa are associated with relevant production traits. during the previous week, plus 10% to account for a feed Many OTU are only present in a small percentage of intake increase as the animal grows. To prevent a pos- the microbiota samples, which implies overdispersion sible association between cage and maternal effects, a due to an excessive number of zero counts that is not maximum of two kits belonging to the same litter were appropriately adjusted with a linear model. Thus, a zero- assigned to the same cage. The dataset used for this study inflated Poisson (ZIP) model might be more suitable to included information on 425 individuals born from 196 estimate the heritability for these traits [25]. In a ZIP litters and housed in 189 cages. The pedigree included model, a given OTU is either not observed (zero counts) information on 9760 individuals, tracing back 20 genera- with probability p or observed with a number of counts tions of ancestors up to animals born in 2005. Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 3 of 15 Sample collection, DNA extraction and sequencing with default parameters. Then, the contigs with a quality At slaughter, cecal samples from each animal were col- score lower than Q19 were discarded, and the remaining lected in a sterile tube, first kept at 4 °C in the laboratory, ones were assigned to samples using the python script and then stored at − 80 °C. Extraction and amplification split_libraries.py with default parameters. The UCHIME of DNA, Illumina library preparation and sequencing algorithm [30] was used to detect and remove the chi- were performed as described in [6]. To facilitate efficient meric sequences generated during PCR. The filtered lysis, 250 mg of each sample were mechanically lysed in contigs were clustered into OTU with a 97% similarity a FastPrep-24TM Homogenizer (MP Biomedicals, LLC, threshold using the script pick_open_reference_otus.py Santa Ana, CA, USA) at a speed of 6 m/s for 60 s. The with default parameters [31] and the Greengenes refer- kit ZR Soil Microbe DNA MiniPrepTM (ZymoResearch, ence database (version gg_13_5_otus) [32]. The obtained Freiburg, Germany) was used to extract whole genomic OTU table was normalized with the cumulative sum scal- DNA. The integrity and purity of the DNA were meas- ing (CSS) method [33]. Finally, the UCLUST consensus ured with a Nanodrop ND-1000 spectrophotometer taxonomy assigner was used to conduct the taxonomic (NanoDrop products; Wilmington, DE, USA) following assignment of representative sequences of each OTU the protocol of Desjardins and Conklin [27]. The F515Y/ (the most abundant sequence of each cluster) by mapping R926 pair of primers (5′-GTG​YCA​GCMGCC​GCG​ the sequences against the Greengenes reference database GTAA-3′, 5′-CCG​YCA​ATTYMTTT​R AG​TTT-3′) [28] gg_13_5_otus. The raw sequence data were deposited in was used to amplify a fragment of the 16S rRNA gene the sequence read archive of NCBI under the BioProject that included the V4-V5 hypervariable regions. An ini- accession number PRJNA524130. The metadata, OTU tial polymerase chain reaction (PCR) was conducted for table, and corresponding taxonomic assignments are each sample with 12.5 µL 2× KAPA HiFi HotStart Ready in Additional file 1: Table S1, Additional file 2: Table S2 Mix, 5 µL forward primer, 5 µL reverse primer and 2.5 and Additional file 3: Table S3, respectively (see Addi- µL template DNA (5 ng/ µL) under the following condi- tional file 1: Table S1, Additional file 2: Table S2 and tions: initial denaturation for 3 min at 95 °C, 25 cycles of Additional file 3: Table S3). These data are also available 30 s at 95 °C, 30 s at 55 °C and 30 s at 72 °C; and final in Qiita (https://​qiita.​ucsd.​edu/) under study ID 14485. extension for 2 min at 72 °C. Then, sequencing adaptors After the bioinformatic processing, 984 representative tiplex ­Nextera® XT kit (Illumina, Inc., San Diego CA, and eight nucleotide dual-indexed barcodes of the mul- traits of the rabbit intestinal microbiota were defined and analyzed in the present study. These microbial traits can USA) were added in a second PCR reaction with 25 µL be categorized into three groups: the relative abundances 2× KAPA HiFi HotStart Ready Mix, 5 µL index i7, 5 µL of 29 genera, 951 CSS-normalized OTU, and four micro- index i5, 10 µL PCR Grade water and 5 µL concentrated bial alpha-diversity indices computed at 10,000 contigs amplicons of the initial PCR. Conditions for this second (total number of OTU observed, Chao1, Shannon and PCR were: an initial denaturation for 3 min at 95 °C, 8 Simpson’s inverse). Relative abundances of genera and cycles of 30 s at 95 °C, 30 s at 55 °C and 30 s at 72 °C; and microbial alpha-diversity indices were standardized by a final extension for 5 min at 72 °C. The libraries obtained subtracting their mean and dividing by their standard were cleaned up with AMPure XP beads, validated by deviation. Finally, these standardized microbial traits and running 1 µL of a 1:50 dilution on a Bioanalyzer DNA CSS-normalized OTU were multiplied by 100 and sub- 1000 chip (Agilent Technologies, Inc., Santa Clara, CA, sequently rounded to the nearest integer. This transfor- USA) to verify their size, and quantified by fluorometry mation was necessary to generate integers that can be with the PicoGreen dsDNA quantification kit (Invitro- treated as counts, thus enabling the adjustment of a ZIP gen, Life Technologies, Carlsbad, CA, USA). After size model. verification, libraries were pooled at equimolar concen- trations and paired-end sequenced in five parallel plates Statistical models on an Illumina MiSeq 2× 250 platform at the Genomics Zero‑inflated  Poisson (ZIP)  mixed model and Bioinformatics Service of the Autonomous Univer- Let y = y1 , y2 , . . . , yn ′ be the vector of the records of sity of Barcelona. a specific microbial trait on n individuals. Zero-infla- tion occurs with probability p, and data for animal i follow a Poisson distribution with parameter i occur Bioinformatics processing of microbial traits with probability (1 − p). Thus, the  probability of observ- The pipeline of the QIIME software (version 1.9.0) [29] ing a zero count is p yi = 0 = p + (1 − p)e−i , and that was used for sequence processing is fully described the probability of observing a count of k (k = 1, 2, …, in [6]. Briefly, paired-end reads were assembled into con- ∞) is p yi = k =   (1−p)e−i ik . Therefore, p is a popula- k! tigs using the python script multiple_join_paired_ends.py tion parameter while i is an animal-specific parameter. Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 4 of 15 Conditioning on both p and , which is the vector of the constraint implies that the density of the resulting distri- i of all the individuals, the likelihood function can be bution will be just 1/6 of the 3-variate uniform distribu- expressed as follows: tion without constraint. As explained later, the marginal   prior distributions are needed to define the BF for the    (1 − p)e−i i yi ratios of variance. Assuming prior independence, the p + (1 − p)e−i  p y|, p = . yi =0 yi >0 yi ! marginal priors of h2, l2 and c2 from the constrained joint (1) distribution are the following Beta distributions: Considering these two re-parameterizations:    2 p h2 = Beta(1, 3) = 3 h2 − 6h2 + 3, (8) i∗ = log(i ), (2)    2 p l2 = Beta(1, 3) = 3 l2 − 6l2 + 3, (9)   p p∗ = log , (3) 1−p  2 the previous conditional likelihood can be expressed as:   p c2 = Beta(1, 3) = 3 c2 − 6c2 + 3. (10)     ∗ ∗  1  ∗ p −exp ( ∗ )  A uniform distribution along the positive real numbers p y| , p =  ∗  e +e i yi =0 1 + ep was also assumed for σP2 . ZA , ZL and ZC are design matri-   ∗ ∗  ces relating the observations with animals in the pedigree  1 e−exp(i )+i yi (9760 levels), litters (196 levels) and cages (189 levels),  ∗  , 1 + ep yi >0 yi ! respectively, and A is the numerator relationship matrix (4) [34]. Uniform priors were assumed for the elements of β exp(p∗ ) and p∗, bounded between − 5 and + 5 for the latter. since i = exp and p = i∗ 1+exp(p∗ ).   The posterior density can be written as: At a subsequent hierarchical level, different factors can be included as a linear model to explain the vector ∗, p(∗ , p∗ , V, β|y) ∝ p(y|∗ , p∗ )p ∗ |V, β p p∗ p(V)p(β)     thus, the assumed distribution for ∗ is the following nor- mal density: p ∗ , p∗ , V, β|y   p ∗ |V, β ∼ MVN(Xβ, V),   (5)     1  ∗ p −exp ( ∗ )  ∝  ∗  e +e i where β is a vector of systematic factors including the yi =0 1 + ep effects of the different categories of the combination   ∗ ∗  1 e−exp(i )+i yi between breeding farm, diet, and feeding regime (6 lev-  ×  ∗  els), of batch (5 levels) and of body weight at weaning (2 yi >0 1 + ep yi ! levels). X is a design matrix that relates the observations  1  × |V|n/2 exp − ∗ − Xβ ′ (V)−1 ∗ − Xβ    to the systematic effects, and V is the covariance matrix 2 between the elements of ∗. The structure of V is not (11) diagonal, and is defined as follows: This model specification is similar to that previously ′ ′ V = σP2 [ZA AZA h2 + ZL ZL l2 + ZC ZC c2 ′ proposed for studying mastitis cases in dairy cows [35].   (6) The differences introduced in our study refer to the + I 1 − h2 − l2 − c2 ], specifications for ∗. In our case, we assume a model in which a number of factors was absorbed into the resid- where σP2 is the phenotypic variance and the scalars h2, l2 ual, while these factors were explicitly fitted in the model and c2 represent the ratios of additive genetic, litter and of Rodrigues-Motta et al. [35] as part of the vector of the cage variances to σP2 . The assumed joint prior distribution means. Apart from these differences in the prior assump- of these ratios was uniform, with the constraint that their tion, the two models are equivalent. The parameteriza- sum must be less than 1: tion used here allows the computation of the BF for the   ratio of variances in a parametric space defined between p h2 , l2 , c2 = 6, if h2 + l2 + c2 ∈ [0, 1], or 0 otherwise. 0 and 1, including both limits, although it is more com- (7) putationally demanding than that of Rodrigues-Motta To guarantee a density of 1, the joint prior distribu- et al. [35]. The reasons are that, first, V must be updated tion under such a constraint was set to 6. Note that the and inverted repeatedly and second, Metropolis-Hasting Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 5 of 15 steps are needed to update the conditional posterior dis- hypotheses were independently tested by computing the tribution of the ratios. In contrast to the case where the BF of h2 = 0 against h2 = 0 (BFh2), l2 = 0 against l2 = 0 effects are explicitly considered in the model [35], the BF (BFl2), and c2 = 0 against c2 = 0 (BFc2). can be computed to test whether the additive genetic, 3 3 litter, and cage effects are null or not, since this model BFh2 =  2  ; BFl2 =   ; and parameterization allows a null value of the ratio. The der- p h = 0|y p l2 = 0|y ivation of the conditional posterior distributions can be 3 BFc2 =  2 followed in the studies of Rodrigues-Motta et al. [35] and  p c = 0|y Varona et al. [36] on which our model is based on. (13) The derivations of these definitions of the BF are reported Linear mixed model (LMM) in Varona et al. [36]. They show how the BF for nested This model can be considered as a simplification of the models is the ratio between the marginal prior density eval- previous one since the generation process assumed for all uated at the value of interest, zero in our case, divided by the data was the same as that assumed for the logarithm the marginal posterior density at the  value of interest (i.e., of the vector of  parameters of the individual Poisson p h2 =0 for the heritability: BFh2 = p h2 =0|y ). The marginal prior distributions (∗) corresponding to the records with non- ( ) zero counts for each trait (transformed CSS-normalized densities under our constrained joint prior distribution of OTU counts, transformed relative abundances of genera, the ratios are Beta(1,3) distributions, which are equal to 3 and transformed alpha-diversity indices). Thus, the dis- at the 0 value for the ratio of variance (Eq. (8)). The evalua- tribution of the data given the model parameters can be tion of the marginal posterior of the ratios at 0 implies that, written as: since these marginal posteriors are only defined up to pro- portionality, the computation of this proportionality is con-    h2 =1  p y|V, β ∼ MVN(Xβ, V). (12) stant: h2 =0 p h2 |y dh2. This integral can be solved  numerically at each iteration. The different BF can be com- Notation and model specifications, including both the puted as follows from the Markov chain Monte Carlo structure of V (Eq. (6)) and the prior definitions (Eq. (7)), (MCMC) output: were assumed to be the same as in the ZIP model, except for the parameters that are specific to the ZIP model. For 3 BFh2 =  , implementation, we used the conditional posterior distri- N p(h2 =0|y )j (14) butions of this LMM derived by Varona et al. [36] since j=1 N we assumed the same prior distributions as theirs. where N is the number of MCMC iterations and p h2 = 0|y j is the evaluation of the marginal posterior  Criteria for comparison of models Two criteria for model choice were applied for each of density of h2 at zero at each iteration j of the sampling the 984 microbial traits analyzed in this study. First, we procedure, which is computed as stated above: evaluated whether the trait was better adjusted with the p h2 = 0|y j   LMM or the ZIP model with the preferred model being  h2 =1  . (15) that with the lowest deviance information criterion 2 |y dh2  2 h =0 p h j (DIC) value [37]. The statistical relevance of the addi- tive genetic, litter and cage effects was evaluated for All the operations were done on the logarithmic scale each trait with the model that best fitted the trait (either and after having saved the evaluations of the marginal pos- LMM or ZIP) using the BF. Thus, for each trait, three BF terior at zero along the MCMC chain to avoid numerical were computed to assess the null hypotheses of whether instabilities during their computation. In this way, it was the additive genetic (a), litter (l) or cage (c) effects have a possible to adjust the evaluations of the marginal poste- null effect versus the alternative hypothesis that assumed rior at zero for their maximum, thus reducing the needed that these factors have a non-null effect. These three numerical accuracy:  � � � � �� �  N p h2 = 0|y � � N 2 = 0|y  j=1 exp log p h j − m  j (16) � = exp log +m , N  N  j=1 Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 6 of 15 where m is the maximum value of the  vector composed Table 2 Phenotypic summary of genera and alpha diversity of the N evaluations of p h2 = 0|y on the logarithmic indices scale. See Sorensen and Gianola [38] for further details. Trait Mean SD BF values were classified according to four levels of evi- dence [39]: BF < 3.2: the denominator model is supported; Genus Methanobrevibacter, % 0.19 0.23 3.2 ≤ BF < 10: substantial evidence favoring the numera- Genus Adlercreutzia, % 0.95 0.43 tor model; 10 ≤ BF < 100: strong evidence favoring the Genus Bacteroides, % 1.65 0.76 numerator model; and BF ≥ 100: decisive evidence favor- Genus Parabacteroides, % 0.21 0.18 ing the numerator model. The reciprocal of the BF value Genus Rikenella, % 0.35 0.24 (1/BF) was used to assess the evidence for the null hypoth- Genus Butyricimonas, % 0.20 0.19 esis, which assumed that the additive genetic, litter or cage Genus Odoribacter, % 0.27 0.22 effects on the trait were null. Genus Clostridium, % 1.05 0.26 Genus Dehalobacterium, % 0.08 0.09 MCMC Bayesian implementation Genus Anaerofustis, % 0.11 0.07 MCMC Bayesian procedures were used to obtain sam- Genus Anaerostipes, % 0.16 0.08 ples from the marginal posterior distributions. This algo- Genus Blautia, % 2.94 0.65 rithm was implemented in a Fortran 90 software, which Genus Butyrivibrio, % 0.11 0.07 is available in our GitHub repository (https://​github.​com/​ Genus Coprococcus, % 2.02 0.42 juanp​ablo-​sanch​ez/​BF-​ZIP). For both the LMM and ZIP Genus Dorea, % 0.47 0.12 models, chains of 10,000 samples were run discarding the Genus Epulopiscium, % 0.11 0.11 first 1000 to allow the algorithm to reach convergence Genus Ruminococcus, % 0.16 0.07 to the marginal posterior distributions. Convergence Genus rc4-4, % 0.15 0.07 diagnostics of the Markov chains was performed by the Genus Faecalibacterium, % 0.20 0.05 Geweke test function with the coda R package [40]. The Genus Oscillospira, % 2.26 0.58 z-statistics of this test for the genetic, litter and cage vari- Genus Phascolarctobacterium, % 0.21 0.24 ance ratios computed for the 984 microbial traits fitted Genus Coprobacillus, % 0.19 0.24 with the LMM or ZIP models (according to DIC) are in Genus p-75-a5, % 0.10 0.07 Additional file 4: Fig. S1 and Additional file 5: Fig. S2, Genus Oxalobacter, % 0.11 0.06 respectively. Although the computational demand of the Genus Desulfovibrio, % 0.46 0.31 parameterization on the variance ratios is high, it allows Genus Campylobacter, % 0.07 0.08 a good mixing. Thus, a reduced number of iterations was Genus Ruminococcus, % 4.32 0.85 needed to properly reach convergence and characterize Genus Anaeroplasma, % 0.20 0.17 the marginal posterior distributions. Genus Akkermansia, % 1.47 0.50 Number of OTU observed 551.05 91.94 Results Shannon 5.07 0.30 Cecal microbial composition and diversity Simpson 0.98 0.01 After bioinformatic sequence processing, we identified Simpson’s inverse 71.01 20.20 951 OTU that were present in at least 5% of the animals. SD: standard deviation Each sample had on average 677 OTU (range from 197 to 841). Table 1 shows the OTU frequencies across rabbit samples. Table 1 OTU frequencies across rabbit cecal samples The taxonomic assignment of representative sequences Frequency (%) Number of these OTU against the Greengenes reference database of OTU gg_13_5_otus (see Additional file 3: Table S3) revealed ≥ 5 to ≤ 10 13 the presence of 29 known genera. Among these, four > 10 to ≤ 25 43 were present in 50 to 75% of the rabbit samples and 25 in > 25 to ≤ 50 121 a minimum of 75% of the animals. Table 2 shows a phe- > 50 to ≤ 75 277 notypic summary of the relative abundances of the 29 > 75 to ≤ 100 497 genera together with the four microbial alpha-diversity indices computed. Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 7 of 15 Adjustment of microbial traits with the LMM versus ZIP variance ratios for the OTU that, according to DIC, are model better adjusted with the LMM, and for those for which According to DIC, the adjustment of all relative abun- the ZIP model is preferable are shown in Figs. 1 and 2, dances of genera and microbial alpha-diversity indices respectively. The same plots corresponding to relative was better with the LMM (lower DIC values) than with abundances of genera and microbial alpha-diversity the ZIP model (see Additional file 6: Table S4). Among indices are shown in Figs. 3, and 4, respectively. In all of the 951 CSS-normalized OTU also analyzed in our the plots, microbial traits are categorized by their fre- study, those with a frequency higher than 25% were bet- quency among samples and by the BF that represent the ter adjusted with the LMM and those with a frequency levels of evidence in favor of the model that included lower than 10% were better adjusted with the ZIP model. additive genetic (a), litter (b), or cage (c) effects. Over- Among the 43 OTU with a frequency between 10 and all, these figures show that the BF did not provide evi- 25%, 18 and 25 OTU were better adjusted with the LMM dence of genetic, litter, or cage effects for ~ 80% of the and the ZIP model, respectively (see Additional file 6: microbial traits analyzed. However, for the traits that Table S4). were declared to be affected by the host genetics, lit- ter or cage (see Additional file 7: Tables S5, S6 and S7, Effects of the host genetics, litter and cage on the rabbit respectively), the magnitude of the estimates of vari- cecal microbiota ance ratios was moderate to high with minimum values Box and whisker plots of the estimated marginal pos- of 0.15 to 0.20. terior means of the heritability, and litter and cage Fig. 1 Marginal posterior of heritability (a), litter (b) and cage (c) variance ratios for OTU adjusted with the normal linear mixed model. OTU are categorized by their frequency among the rabbit samples. Each OTU is represented by a dot colored in red, yellow, blue or green based on Bayes factor for no, substantial, strong or decisive evidence, respectively, of additive genetic (a), litter (b) or cage (c) effects Fig. 2 Marginal posterior of heritability (a), litter (b) and cage (c) variance ratios for OTU adjusted with the zero-inflated Poisson model. OTU are categorized by their frequency among the rabbit samples. Each OTU is represented by a dot colored in red, yellow, blue or green based on Bayes factor for no, substantial, strong or decisive evidence, respectively, of additive genetic (a), litter (b) or cage (c) effects Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 8 of 15 Fig. 3 Marginal posterior of heritability (a), litter (b) and cage (c) variance ratios for relative abundances of genera adjusted with the normal linear mixed model. Genera are categorized by their frequency among the rabbit samples. Each genus is represented by a dot colored in red, yellow, blue or green based on Bayes factor for no, substantial, strong or decisive evidence, respectively, of additive genetic (a), litter (b) or cage (c) effects Fig. 4 Marginal posterior of heritability (a), litter (b) and cage (c) variance ratios for microbial alpha-diversity indices adjusted with the normal linear mixed model. Each index is represented by a dot colored in red, yellow, blue or green based on Bayes factor for no, substantial, strong or decisive evidence, respectively, of additive genetic (a), litter (b) or cage (c) effects Table 3 Mean (standard deviation) of Bayes factors and heritability estimates for OTU under genetic control adjusted with the normal LMM Frequency (%) Substantial evidence of genetic Strong evidence of genetic influence (10 Strong evidence of genetic influence (3.2 ≤ BFh2 < 10) ≤ BFh2 < 100) influence (BFh2 ≥ 100) BFh2 h2 n BFh2 h2 n BFh2 h2 n All 5.31 (1.62) 0.21 (0.03) 108 24.83 (18.09) 0.26 (0.05) 33 159.38 (15.27) 0.19 (0.04) 2 > 10 to ≤ 25 5.03 (1.32) 0.20 (0.03) 6 26.92 (12.50) 0.23 (0.03) 4 170.18 (–) 0.16 (–) 1 > 25 to ≤ 50 5.49 (1.79) 0.21 (0.03) 23 36.19 (31.80) 0.25 (0.04) 6 – – 0 > 50 to ≤ 75 5.46 (1.70) 0.21 (0.04) 44 21.23 (12.57) 0.28 (0.05) 12 148.58 (–) 0.22 (–) 1 > 75 to ≤ 100 5.07 (1.48) 0.19 (0.03) 35 21.82 (14.73) 0.27 (0.07) 11 – – 0 BFh2 = Bayes factor of the model with additive genetic effects against the same model without additive genetic effects Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 9 of 15 The results summarized in Figs. 1, 2, 3 and 4 are pre- Table 5 Bayes factors, marginal posterior means (standard sented in the following paragraphs of the “Results” deviations) of the heritability estimates for genera and alpha- section. Bayes factors, marginal posterior means and diversity indices influenced by genetic effects standard deviations of variance ratios for the OTU that, Trait ­2 Mean (SD) h P ­(h2 > 0.1) BFh2 based on the BF, were declared to be influenced by host genetics, litter, or cage effects (together with the associ- Genus Methanobrevibacter 0.21 (0.13) 0.79 7.75 ated probability of these estimates being greater than Genus Butyricimonas 0.27 (0.19) 0.79 4.39 0.10) are included in Additional file 8: Table S8 (heritable Genus Odoribacter 0.19 (0.13) 0.71 3.50 OTU), Additional file 9: Table S9 (OTU influenced by the Genus Bacteroides 0.29 (0.17) 0.87 13.88 nursing environment) and Additional file 10: Table S10 Genus Parabacteroides 0.35 (0.17) 0.91 31.15 (OTU influenced by cage effects). Genus Dehalobacterium 0.29 (0.19) 0.83 8.62 Genus Blautia 0.20 (0.12) 0.78 7.01 Microbial traits under genetic control Genus Epulopiscium 0.17 (0.11) 0.70 5.85 Table 3 summarizes the marginal posterior means of Genus Phascolarctobacterium 0.19 (0.12) 0.73 10.22 the heritability estimates for the OTU, categorized by Genus Oxalobacter 0.21 (0.13) 0.78 6.12 frequency, which were better adjusted with the normal Number of OTU observed 0.28 (0.17) 0.84 7.30 LMM according to DIC and for which the BF provided Shannon 0.18 (0.13) 0.70 3.41 evidence in favor of host genetics control. Similarly, SD: standard deviation; BFh2: Bayes factor of the model with additive genetic effects against the same model without additive genetic effects Table 4 summarizes the marginal posterior means of the heritability estimates for OTU, categorized by frequency, which were better adjusted with the ZIP model according to DIC and for which the BF provided evidence in favor noted that these estimates have large standard deviations of host genetics control. as a consequence of our limited sample size. Neverthe- The BF provided some level of evidence in favor of less, it is worth mentioning that for 51 of the 154 OTU host genetics control for 154 of the 951 OTU analyzed. identified as being under genetic control, the probability The BF between models with and without an additive that their heritability estimate is higher than 0.10 was genetic effect provided evidence of a substantial (BFh2 ≥ equal or greater than 0.80. 3.2) genetic control for 108 and 10 of the OTU that were Table 5 shows the marginal posterior means and better adjusted with the normal LMM and ZIP mod- standard deviations of the heritability estimates, els, respectively. Evidence of a strong (10 ≤ BFh2 < 100) together with the probability of these estimates being genetic control was found for 33 and one OTU that were greater than 0.10, for relative abundances of genera better adjusted with the normal LMM and ZIP models, and alpha-diversity indices. The BF provided evidence respectively. Finally, decisive (BFh2 ≥ 100) evidence of in favor of genetic control for the relative abundances genetic control was found for two OTU that were better of 10 genera. This evidence was substantial ( BFh2 ≥ 3.2) adjusted with the normal LMM. The taxonomic assign- for the genera Dehalobacterium, Epulopiscium, Metha- ment of these two OTU revealed that one belongs to the nobrevibacter, Butyricimonas, Odoribacter, Blautia and genus Bacteroides and the other to the genus Parabacte- Oxalobacter; and strong (10 ≤ BFh2 < 100) for the gen- roides, and the marginal posterior means (standard devi- era Phascolarctobacteirum, Bacteroides and Parabacte- ations) of their heritability estimates were 0.16 (0.07) and roides. The estimates (marginal posterior means) of the 0.22 (0.08), respectively (see Additional file 8: Table S8). heritability for these genera ranged from 0.17 to 0.35. Overall, the estimates of the heritability for these OTU The highest heritability estimates, accompanied by were moderate (from 0.12 to 0.40). However, it should be high BF values, were found for the genera Bacteroides, Table 4 Mean (standard deviation) of Bayes factors and heritability estimates for OTU under genetic control adjusted with the ZIP model Frequency (%) Substantial evidence of genetic influence (3.2 ≤ BFh2 < 10) Strong evidence of genetic influence (10 ≤ BFh2 < 100) BFh2 h2 n BFh2 h2 n All 4.49 (1.73) 0.24 (0.04) 10 11.37 (–) 0.18 (–) 1 ≥ 5 to ≤ 10 3.62 (0.35) 0.24 (0.01) 5 11.37 (–) 0.18 (–) 1 > 10 to ≤ 25 5.37 (2.16) 0.23 (0.06) 5 – – 0 BFh2= Bayes factor of the model with additive genetic effects against the same model without additive genetic effects Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 10 of 15 Parabacteroides and Dehalobacterium, for which the The BF provided some level of evidence in favor of a lit- probability that their heritability estimates are higher ter effect for 215 of the 951 OTU analyzed. Six of them than 0.10 was greater than 0.80. Although the evidence were better adjusted (lower DIC values) with the ZIP in favor of genetic control was strong for the genus model and the remaining 209 were better adjusted with Phascolarctobacterium, its heritability (0.19) estimate the LMM. BF values between models with and without was not one of the highest (P(h2 > 0.1) = 0.73). However, a litter effect provided evidence of a substantial (BFl2 1/BFh2 values higher than 3.2 were obtained only for 10 ≥ 3.2) litter effect for 81 and three OTU that were bet- OTU and the genus Coprococcus, which are clearly not ter adjusted with the normal LMM and the ZIP model, heritable. respectively. Strong evidence (10 ≤ BFl2 < 100) of a litter Finally, regarding the traits that globally integrate the effect was found for 63 and three OTU that were bet- rabbit cecal microbiota, substantial evidence of genetic ter adjusted with the normal LMM and the ZIP model, control was found for number of OTU observed and respectively. Finally, decisive (BFl2 ≥ 100) evidence of Shannon indices. The highest heritability estimate was a litter effect was found for 65 OTU that were better found for the number of OTU observed index ­(h2 = 0.28; adjusted with the normal LMM. The taxonomic assign- BFh2 = 7.30). The additive genetic background for these ment of these OTU revealed that most of them belong to two traits was clearly demonstrated by a probability the genera Parabacteroides, Phascolarctobacterium, and greater than 0.70 that their heritability estimates are the species eggerthii and fragilis of the genus Bacteroides. higher than 0.10. Overall, the marginal posterior means of the litter vari- ance ratio ranged from 0.12 to 0.19 (Tables 6 and 7) but Microbial traits influenced by the nursing environment the estimates for this ratio reached values ranging from The marginal posterior means of litter variance ratio 0.37 to 0.54 for the OTU for which large BF values were for OTU, categorized by frequency, which were better observed (see Additional file 9: Table S9). Eighty-nine of adjusted with the normal LMM according to DIC and for the 215 OTU that were declared to be influenced by litter which the BF provided evidence in favor of a litter effect effects had a probability equal or greater than 0.80 that are summarized in Table 6, and the corresponding ones their litter variance ratio be greater than 0.10. The mar- which were better adjusted with the ZIP model according ginal posterior means of the litter variance ratio were to DIC are in Table 7. greater than 0.50 for 12 OTU of which six belong to the Table 6 Mean (standard deviation) of Bayes factors and litter variance ratio estimates for OTU influenced by a litter effect adjusted with the normal LMM Frequency (%) Substantial evidence of litter Strong evidence of litter influence (10 Decisive evidence of litter influence influence (3.2 ≤ BFl2 < 10) ≤ BFl2 < 100) (BFl2 ≥ 100) BFl2 l2 n BFl2 l2 n BFl2 l2 n All 5.29 (1.79) 0.12 (0.01) 81 29.46 (22.86) 0.17 (0.02) 63 ∞ (∞) 0.37 (0.15) 65 > 10 to ≤ 25 4.08 (0.80) 0.12 (0.00) 2 10.09 (–) 0.14 (–) 1 ∞ (∞) 0.54 (0.17) 9 > 25 to ≤ 50 5.45(1.73) 0.12 (0.01) 16 25.23 (16.74) 0.17 (0.02) 13 ∞ (∞) 0.40 (0.14) 19 > 50 to ≤ 75 5.49 (1.90) 0.12 (0.01) 28 34.90 (26.65) 0.17 (0.02) 23 ∞ (∞) 0.33 (0.12) 27 > 75 to ≤ 100 5.12 (1.77) 0.12 (0.01) 35 27.51 (21.94) 0.17 (0.02) 26 7.98E4 (2.33E5) 0.27 (0.06) 10 BFl2= Bayes factor of the model with litter effects against the same model without litter effects Table 7 Mean (standard deviation) of Bayes factors and litter variance ratio estimates for OTU influenced by a litter effect adjusted with the ZIP model Frequency (%) Substantial evidence of litter influence (3.2 ≤ BFl2 < 10) Strong evidence of litter influence (10 ≤ BFl2 < 100) 2 BFl2 l N BFl2 l2 n All 3.77 (0.46) 0.19 (0.00) 3 36.69 (23.57) 0.46 (0.14) 3 ≥ 5 to ≤ 10 3.80 (0.64) 0.19 (0.00) 2 59.05 (–) 0.62 (–) 1 > 10 to ≤ 25 3.70 (–) 0.19 (–) 1 25.51 (19.01) 0.38 (0.04) 2 BFl2= Bayes factor of the model with litter effects against the same model without litter effects Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 11 of 15 Table 8 Bayes factors, marginal posterior means (standard Previous studies in humans and different livestock spe- deviations) of litter variance ratio for genera and alpha-diversity cies have pointed out the existence of host genetics con- indices influenced by litter effects trol of the gut microbiota, but, to our knowledge, this has Trait Mean (SD) l­ 2 P(l2 > 0.1) BFl2 not been investigated in the rabbit. For the first time, we have evaluated the host genetics, litter, and cage effects Genus Butyricimonas 0.28 (0.10) 0.96 728.23 on the microbial composition of the rabbit cecum, which Genus Odoribacter 0.14 (0.08) 0.64 7.21 is the organ that contains the greatest microbial diver- Genus Bacteroides 0.27 (0.09) 0.97 809.53 sity and complexity [2]. In this study, we defined a set of Genus Parabacteroides 0.47 (0.10) 1.00 1.50E11 984 microbial traits that represent cecal microbial com- Genus Rikenella 0.32 (0.08) 1.00 3.67E4 position and diversity with different levels of complexity. Genus Dehalobacterium 0.37 (0.10) 1.00 9.66E4 Although genus assignment of every queried sequence Genus Anaerofustis 0.15 (0.08) 0.68 8.80 would be the desired outcome, it should be noted that Genus Epulopiscium 0.12 (0.07) 0.58 4.52 the resolution of the 16S rDNA locus with the MiSeq Genus Phascolarctobacterium 0.66 (0.07) 1.00 ∞ technology used in this study allowed annotation of only Genus Desulfovibrio 0.17 (0.09) 0.78 16.53 21% of all the sequences at the genus level. Nevertheless, Genus Campylobacter 0.16 (0.08) 0.73 11.91 given the importance of the functional roles played by Number of OTU observed 0.17 (0.09) 0.75 15.55 bacteria that can be assigned at the genus level, it is rele- SD: standard deviation; BFl2: Bayes factor of the model with litter effects against vant to provide estimates that inform about the influence the same model without litter effects of genetic, litter and cage effects for groups of reads that belong to the same genus. The CSS-normalized abundances of 951 OTU can be genus Bacteroides, four to the genus Phascolarctobacte- considered as the most specific level of defining a micro- rium, one to the genus Parabacteroides and one to the bial community, but their frequencies vary widely among genus Rikenella (see Additional file 9: Table S9). It should samples. This means that some OTU are present in all or be mentioned that the 1/BFl2 values were greater than 3.2 almost all the animals (i.e., core OTU), while others are for 107 OTU, which are not influenced by a litter effect. only detected in some animals. The distribution of the The marginal posterior means of the litter variance OTU that are present only in a small percentage of the ratio, together with the associated probability of these analyzed animals was clearly far from a normal distribu- estimates being greater than 0.10, for relative abundances tion and, thus not surprisingly, these were better adjusted of genera and alpha-diversity indices are in Table 8. The with the ZIP model. In spite of this, all the microbial number of OTU observed index was found to be influ- traits analyzed in this study with a frequency higher than enced by the nursing environment (­l2 = 0.17; BFl2 = 25% were better adjusted with the normal LMM model 15.55). The BF also provided evidence of a litter effect according to DIC. For the traits with a clear excess of for 10 genera. Very strong evidence of a litter effect was zero counts (i.e., having a frequency among samples found for the genera Butyricimonas ­(l2 = 0.28), Bacte- lower than 15%), DIC clearly favored the ZIP model. Pre- roides ­(l2 = 0.27), Parabacteroides ­(l2 = 0.47), Rikenella vious microbiome studies have also modeled microbiome ­(l2 = 0.32), Dehalobacterium ­(l2 = 0.37) and Phascolarcto- data with the ZIP model to account for the excess of zero bacterium ­(l2 = 0.66) with decisive BF values (BFl2 ≥ 100) counts of many taxa that are rare and only detected in and P(l2 > 0.1) = 0.96. However, the genera Coprococcus, a small proportion of the samples [24, 41]. Such studies rc4-4 and Faecalibacterium are not influenced by a litter argue that the application of a conventional linear model effect (1/BFl2 > 3.2). is inappropriate for zero-inflated data. In our study, the ZIP model performed better than the LMM only for the Discussion microbial traits with a very marked excess of zero counts. The influence of many external factors on the rabbit cecal However, it should be noted that our choice of the ZIP microbial composition and diversity is well documented model, as an alternative to the LMM, was based only on [4–8], but not that of the host genetics. To shed light on its ability to handle the excess of zero counts, and not on this, we report the heritability estimates, together with whether the non-zero inflated component of the model the litter and cage (results for cage effects can be found was properly adjusted to a Poisson distribution. A lack in Additional file 7: Tables S5, S6, and S7) variance ratios of adjustment of transformed CSS records to the Poisson estimates, for microbial traits that were fitted with LMM model could perfectly explain the fact that the ZIP model and ZIP mixed models. Moreover, we have assessed the performed better than the LMM only for the micro- statistical relevance of these estimates through BF. bial traits that have a very large excess of zero counts. Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 12 of 15 Although we have derived some biological conclusions tract of mammals is considered to occur immediately from the application of our models, we are aware of the after birth when the newborns acquire their immature limitation due to the application of LMM on records with microbiota from a combination of maternal and external a marked excess of zero counts. Handling such issues of microbes [3, 48]. The impact of the nursing environment CSS counts remains a challenge. on the relative abundances of these genera still prevails The BF provided evidence in favor of host genetics con- at slaughter age when we collected the cecal samples trol for 34% and 16% of the genera and OTU present in for this study. Remarkably, the ratio of phenotypic vari- the rabbit cecum, respectively. These results are in line ance due to litter effects exceeds 0.50 for six OTU that with the analysis of heritability estimates conducted in belong to the genus Bacteroides and for one OTU that is humans by Goodrich et al. [42] who found evidence of a taxonomically assigned to the genus Parabacteroides. It is genetic control for 10% of their 945 identified taxa, and also worth noting that cage seems to have an important also with an assessment in cattle of the host genetics effect in the relative abundances of some members of the influence on the rumen microbiota [13] that showed that family S24-7. Bacteria within this family, which is part of 34% of the microbial taxa analyzed (from genus to phy- the order Bacteroidales, have been shown to be domi- lum levels) were heritable. Our heritability estimates for nant in the mouse gut microbiota and have been detected the relative abundances of the genera and OTU that were in the gastrointestinal tract of different mammals. The declared to be under host genetics control by the BF were classification of this family remains ambiguous because moderate. This is also in agreement with earlier studies in its members have not been cultured, but the functional humans and other livestock species. However, it is note- analysis conducted by Lagkouvardos et al. [49] renamed worthy that these studies suggested that the main herit- it as the family Muribaculaceae. A recent study on mice able bacteria belong to the phylum Firmicutes, whereas showed that members of the family Muribaculaceae are taxa encompassed by the phylum Bacteroidetes are gen- major mucin monosaccharide foragers in the gut [50]. erally not heritable [13, 43]. A discussion of the results High heritability estimates, accompanied by strong evi- regarding the effects of host genetics, litter, and cage on dence of host genetics control based on the BF, were also taxa that are encompassed by the phyla Bacteroidetes and estimated for the genera Dehalobacterium ­(h2 = 0.29) Firmicutes, as well as on microbial alpha-diversity indi- and Butyricimonas ­(h2 = 0.27). Both genera belong to the ces, will be presented below. phylum Firmicutes and have been previously reported In our study, according to the BF, the strongest evi- as heritable in humans [42, 51]. Such studies reported a dence of host genetics control concerned two OTU that module of co-occurring heritable families within which are taxonomically assigned to the genera Bacteroides and the family Christensenellaceae was the hub (i.e., the node Parabacteroides, which are both part of the phylum Bac- connected to most other nodes) connected to heritable teroidetes. Moreover, the heritability estimates were high- families Methanobacteriaceae and Dehalobacteriaceae. est for these two genera: 0.35 and 0.29 for Parabacteroides Interestingly, we also found substantial evidence of and Bacteroides, respectively. In pigs, Chen et al. [17] and genetic control for the genus Methanobrevibacter, which Bergamaschi et al. [44] also reported some heritable taxa is encompassed by the family Methanobacteriaceae. that belong to the phylum Bacteroidetes. The species that Methanobrevibacter is the single genus of the phylum belong to the genera Bacteroides and Parabacteroides are Euryarchaeota that was detected in the rabbit cecum. It anaerobic Gram-negative bacteria, which are involved in encompasses different hydrogenotrophic methane-pro- the degradation of plant polysaccharides and amino acid ducing species whose abundances have been associated fermentation, amino acid transport, and cell motility in with single-nucleotide polymorphisms that are located the gastrointestinal microbiota of the growing rabbit [45, in a long noncoding RNA, however, this link remains 46]. Although the BF and our heritability estimates for uncertain [52]. In addition, taxa belonging to the family the genera Bacteroides and Parabacteroides clearly reveal Methanobacteriaceae were reported to have heritabil- the existence of a host genetics control, the environmen- ity estimates higher than 0.50 in a beef cattle population tal effect of litter has a profound impact on the relative [53]. It is worth noting that we also found heritability abundances of both genera (­l2 Parabacteroides = 0.47; ­l2 estimates statistically higher than zero for the genera Bacteroides = 0.27). The nursing environment provided Blautia and Odoribacter, which is consistent with previ- by the mother and the siblings also has an important ous results in humans [54, 55]. impact on the relative abundance of the genus Rikenella Our results revealed an important litter effect on ­(l2 Rikenella = 0.32), which is also part of the phylum Bac- the relative abundances of the genera Butyricimonas teroidetes. Litter effects play an important role on pheno- ­(l2 = 0.28), Dehalobacterium ­(l2 = 0.37) and Phascolarcto- typic traits related to rabbit growth and feed efficiency bacterium ­(l2 = 0.66). Estimates of the heritability and BF [47]. The microbial colonization of the gastrointestinal also suggested a genetic control for these three genera, Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 13 of 15 but the effect exerted by the nursing environment seems to reduce this genus could be recommended, it would to be stronger. Regarding the role played by cage environ- only account for a certain part of the methane emissions. mental effects, they were strong for some species encom- Moreover, members of this genus could be beneficial for passed by the genus Ruminococcus. On the contrary, other relevant traits, thus selecting for a given trait via the host genetics and litter effects do not seem to have any microbiota might result in negative responses for other relevant influence on this genus. However, La Reau et al. traits of interest. A direct selection somehow guarantees [56] and Li et al. [13] found that the abundance of the a balanced modification of all the elements involved in genus Ruminococcus was influenced by the host genet- the metabolic pathway of the trait. ics. Remarkably, this genus displays a large diversity and, Finally, it must be noted that the mechanisms underly- our results identified four OTU that are taxonomically ing the host genetics control on cecal microbiota remain assigned as Ruminococcus and are clearly influenced by unknown. Future genome-wide association studies with the host genetics. large datasets are necessary to identify the host genomic On a global scale, our results suggest that a substan- regions that are involved in the control of the overall tial part of the cecal microbial variability is under host microbial diversity and abundances of specific taxa. genetics control since the BF provided evidence in favor of genetic control for the number of OTU observed and Conclusions Shannon indices. In line with previous heritability assess- The Bayesian analysis of a set of 984 microbial traits con- ments of alpha-diversity in pigs [18, 44] and humans [43], ducted in this study with LMM and ZIP mixed models we have found a clear genetic background for the num- has allowed us to dissect the additive genetic, litter and ber of OTU observed index (­h2 = 0.28). Microbial com- cage effects on different levels of complexity of the rab- plexity can be summarized into alpha-diversity indices, bit cecal microbiota based on BF. Fitting these micro- which are heritable traits that could potentially be incor- bial traits with a LMM model was preferable except for porated into breeding programs if an association with a the analyses of the CSS-normalized abundances of rare production trait of economic interest was demonstrated. OTU that are characterized by a marked excess of zero Nonetheless, it is important to bear in mind that alpha- counts and were better adjusted with the ZIP model. The diversity at weaning might not be an accurate predictor calculation of BF as an assessment tool of the statistical of diversity at later stages in rabbit life. Rabbit cecum relevance of the estimates of heritability, litter and cage hosts a rich and complex microbial ecosystem that is variance ratios has allowed us to unravel the different lev- shaped by many non-genetic factors, however, a signifi- els of evidence in favor of these effects on the global cecal cant proportion of the microbial traits analyzed in this microbial composition and on an important proportion study showed moderate heritabilities. Although the cecal of OTU and relative abundances of genera. We found a samples analyzed here were collected from nearly adult strong influence of the host genetics and the nursing rabbits, these estimates should be interpreted with cau- environment for members of the genera Bacteroides and tion since microbial composition changes as the animal Parabacteroides, and of cage effects for the family S24- grows and does not stabilize until the animal reaches 7 and the genus Ruminococcus. Our findings provide adulthood. As we have stated, recent studies in differ- support that the host genetics, cage and nursing envi- ent livestock species have attempted to dissect the host ronment contribute to the variation of the rabbit cecal genetics control of gut microbiota but without paying microbial composition, but functional and genome-wide enough attention to non-genetic factors, such as litter or association studies are needed to decipher the underlying cage effects, which can be even more relevant than the mechanisms. additive genetic effects. For the first time, we have evaluated the role played by Supplementary Information host genetics, litter and cage effects on a set of traits that The online version contains supplementary material available at https://​doi.​ attempt to represent the rabbit cecal microbiota at differ- org/​10.​1186/​s12711-​022-​00738-2. ent levels of complexity. Understanding how these effects Additional file 1: Table S1. Metadata associated with the 425 rabbit cecal influence the intestinal microbiota is relevant from a bio- samples analyzed in this study. logical perspective. One example is the genus Methano- Additional file 2: Table S2. Table of filtered and CSS-normalized OTU. brevibacter, which is clearly heritable and seems to be Additional file 3: Table S3. Taxonomic assignment of representative linked to methane emissions. The host genetics control sequences of each OTU included in Additional file 2: Table S2. for methane emissions and the relative abundance of this Additional file 4: Figure S1. Geweke z-statistics for genetic, litter and genus would offer the possibility to alter the microbial cage variance ratios computed for the microbial traits that were better composition through selection and to breed for rabbits adjusted with the normal linear mixed model. with a reduced impact on climate. Although selection Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 14 of 15 assignments are all included in Additional file 1: Table S1, Additional file 2: Additional file 5: Figure S2. Geweke z-statistics for genetic, litter and Table S2 and Additional file 3: Table S3, respectively. These data are available in cage variance ratios computed for the microbial traits that were better Qiita (https://​qiita.​ucsd.​edu/) [57] under study ID 14485. adjusted with the zero-inflated Poisson model. Additional file 6: Table S4. Deviance information criterion (DIC) values computed for each microbial trait with the normal linear mixed model Declarations and the zero-inflated Poisson model. Ethics approval and consent to participate Additional file 7: Table S5. Mean (standard deviation) of Bayes factors This study was carried out in accordance with the recommendations of the and cage variance ratio estimates for OTU that are influenced by cage animal care and use committee of the Institute for Food and Agriculture effects and adjusted with the normal LMM. Table S6. Mean (standard Research and Technology (IRTA). The protocol was approved by the commit‑ deviation) of Bayes factor and cage variance ratio estimates for OTU that tee of the IRTA. are influenced by cage effect and adjusted with the ZIP model. Tables S7. Bayes factors, marginal posterior means (standard deviations) of cage Consent for publication variance ratio for genera and alpha-diversity indices that are influenced by Not applicable. cage effects. Additional file 8: Table S8. Bayes factors, marginal posterior means and Competing interests standard deviations of heritability estimates for OTU under genetic control The authors declare that they have no competing interests. together with the associated probability of these estimates being higher than 0.10. The frequency among the rabbit samples and the taxonomic Author details 1 assignment are also provided for each OTU. Institute of Agrifood Research and Technology (IRTA)-Animal Breeding and Genetics, Caldes de Montbui, Barcelona, Spain. 2 Centre for Research Additional file 9: Table S9. Bayes factors, marginal posterior means and in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, Bellaterra, Barcelona, standard deviations of litter variance ratio for OTU influenced by litter Spain. 3 Veterinary Faculty, University of Zaragoza, Saragossa, Spain. effects together with the associated probability of these estimates being higher than 0.10. The frequency among the rabbit samples and the taxo‑ Received: 30 September 2021 Accepted: 17 June 2022 nomic assignment are also provided for each OTU. Additional file 10: Table S10. Bayes factors, marginal posterior means and standard deviations of cage variance ratio for OTU influenced by cage effects together with the associated probability of these estimates being higher than 0.10. The frequency among the rabbit samples and the References taxonomic assignment are also provided for each OTU. 1. Flint HJ, Scott KP, Louis P, Duncan SH. The role of the gut microbiota in nutrition and health. Nat Rev Gastroenterol Hepatol. 2012;9:577–89. 2. Gouet PH, Fonty G. Changes in the digestive microflora of holoxenic∗ Acknowledgements rabbits from birth until adulthood. Ann Biol Anim Biochim Biophys. Project funded by Ministry of Science and Innovation-State Research Agency 1979;19:553–66. (AEI, Spain, 10.13039/501100011033) grant number PID2019-108829RB-I00 3. Combes S, Michelland RJ, Monteils V, Cauquil L, Soulié V, Tran NU, et al. to MPE. This work was also supported by grant CEX2019-000902-S funded by Postnatal development of the rabbit caecal microbiota composition and MCIN/AEI/10.13039/501100011033 and by the CERCA Programme/Generalitat activity. FEMS Microbiol Ecol. 2011;77:680–9. de Catalunya (Spain). We would like to thank Oscar Perucho, Josep Ramon 4. Zhu Y, Sun Y, Wang C, Li F. Impact of dietary fibre: starch ratio in shap‑ and Carmen Requena (staff of Unitat de Cunicultura, IRTA) for animal care and ing caecal archaea revealed in rabbits. J Anim Physiol Anim Nutr. their contribution to data recording. We acknowledge Oriol Rafel, Marc Viñas, 2017;101:635–40. Miriam Guivernau and Olga González for their help to collect and process the 5. Chen SY, Deng F, Jia X, Liu H, Zhang GW, Lai SJ. Gut microbiota profiling biological samples. We also acknowledge Armand Sánchez, Nicolas Boulanger with differential tolerance against the reduced dietary fibre level in rabbit. and Joana Ribes (Genomics and NGS Unit, CRAG) for their assistance in mas‑ Sci Rep. 2019;9:288. sive library preparations. 6. Velasco-Galilea M, Guivernau M, Piles M, Viñas M, Rafel O, Sánchez A, et al. Breeding farm, level of feeding and presence of antibiotics in the feed Author contributions influence rabbit cecal microbiota. Anim Microbiome. 2020;2:40. JPS and MP conceived the experimental design. MVG, JPS and MP collected 7. Abecia L, Fondevila M, Balcells J, Lobley GE, McEwan NR. The effect of the biological samples. MVG and MP processed the samples in the laboratory. medicated diets and level of feeding on caecal microbiota of lactating MVG processed and analyzed the sequencing data, interpreted data, and rabbit does. J Appl Microbiol. 2007;103:787–93. prepared the figures and tables. JPS and LV contributed to the analysis of the 8. Zou F, Zeng D, Wen B, Sun H, Zhou Y, Yang M, et al. Illumina Miseq sequencing data. MVG and JPS wrote the manuscript. JPS, MP, YRC and LV platform analysis caecum bacterial communities of rex rabbits fed with helped interpreting the data, and revised the manuscript. All authors read and different antibiotics. AMB Express. 2016;6:100. approved the final manuscript. 9. Goodrich JK, Davenport ER, Clark AG, Ley RE. The relationship between the human genome and microbiome comes into view. Ann Rev Genet. Funding 2017;51:413–33. The experimental design of this work was conducted thanks to funding from 10. Cahana I, Iraqi FA. Impact of host genetics on gut microbiome: take- the INIA project RTA2011-00064-00-00. This study was part of the Feed-a-Gene home lessons from human and mouse studies. Anim Model Exp Med. project that received funding from the European Union’s H2020 program 2020;3:229–36. under grant agreement no. 633531, and the Spanish project RTI2018- 11. Difford GF, Plichta DR, Løvendahl P, Lassen J, Noel SJ, Højberg O, et al. 097610R-I00. MVG was a recipient of a “Formación de Personal Investigador Host genetics and the rumen microbiome jointly associate with methane (FPI)” pre-doctoral fellowship from INIA, associated with the research project emissions in dairy cows. PLoS Genet. 2018;14:e1007580. RTA2014-00015-C2-01. YRC is recipient of a Ramon y Cajal post-doctoral 12. Wallace RJ, Sasson G, Garnsworthy PC, Tapio I, Gregson E, Bani P, et al. fellowship (RYC2019-027244-I) from the Spanish Ministry of Science and A heritable subset of the core rumen microbiome dictates dairy cow Innovation. productivity and emissions. Sci Adv. 2019;5:eaav8391. 13. Li F, Li C, Chen Y, Liu J, Zhang C, Irving B, et al. Host genetics influence the Availability of data and materials rumen microbiota and heritable rumen microbial features associate with The raw sequence data were deposited in the sequence read archive of NCBI feed efficiency in cattle. Microbiome. 2019;7:92. under the accession number SRP186982 (BioProject PRJNA524130). Metadata, 14. Saborío-Montero A, Gutiérrez-Rivas M, García-Rodríguez A, Atxaerandio the filtered and CSS-normalized OTU table and corresponding taxonomic R, Goiri I, Lopez de Maturana E, et al. Structural equation models to Velasco‑Galilea et al. Genetics Selection Evolution (2022) 54:46 Page 15 of 15 disentangle the biological relationship between microbiota and complex 37. Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A. Bayesian measures traits: methane production in dairy cattle as a case of study. J Anim Breed of model complexity and fit. J R Stat Soc B. 2002;64:583–639. Genet. 2020;137:36–48. 38. Sorensen D, Gianola D. Likelihood, Bayesian and MCMC methods in 15. Zhang Q, Difford G, Sahana G, Løvendahl P, Lassen J, Lund MS, et al. quantitative genetics. New York: Springer; 2002. Bayesian modeling reveals host genetics associated with rumen 39. Jeffreys H. The theory of probability. Nature. 1922;109(2727):132–3. microbiota jointly influence methane emission in dairy cows. ISME J. 40. Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and 2020;14:2019–33. output analysis for MCMC. R News. 2006;6:7–11. 16. Camarinha-Silva A, Maushammer M, Wellmann R, Vital M, Preuss S, Benne‑ 41. Lee KH, Coull BA, Moscicki AB, Paster BJ, Starr JR. Bayesian variable selec‑ witz J. Host genome influence on gut microbial composition and micro‑ tion for multivariate zero-inflated models: application to microbiome bial prediction of complex traits in pigs. Genetics. 2017;206:1637–44. count data. Biostatistics. 2020;21:499–517. 17. Chen C, Huang X, Fang S, Yang H, He M, Zhao Y, et al. Contribution of host 42. Goodrich JK, Waters JL, Poole AC, Sutter JL, Koren O, Blekhman R, et al. genetics to the variation of microbial composition of cecum lumen and Human genetics shape the gut microbiome. Cell. 2014;159:789–99. feces in pigs. Front Microbiol. 2018;9:2626. 43. Goodrich JK, Davenport ER, Waters JL, Clark AG, Ley RE. Cross-species 18. Lu D, Tiezzi F, Schillebeeckx C, McNulty NP, Schwab C, Shull C, et al. Host comparisons of host genetic associations with the microbiome. Science. contributes to longitudinal diversity of fecal microbiota in swine selected 2016;352:532–5. for lean growth. Microbiome. 2018;6:4. 44. Bergamaschi M, Maltecca C, Schillebeeckx C, McNulty NP, Schwab C, 19. Ramayo-Caldas Y, Prenafeta-Boldú F, Zingaretti LM, Gonzalez-Rodriguez Shull C, et al. Heritability and genome-wide association of swine gut O, Dalmau A, Quintanilla R, et al. Gut eukaryotic communities in pigs: microbiome features with growth and fatness parameters. Sci Rep. diversity, composition and host genetics contribution. Anim Microbiome. 2020;10:10134. 2020;2:18. 45. Dai ZL, Wu G, Zhu WY. Amino acid metabolism in intestinal bacteria: links 20. Campbell JH, Foster CM, Vishnivetskaya T, Campbell AG, Yang ZK, Wymore between gut ecology and host health. Front Biosci. 2011;16:1768–86. A, et al. Host genetic and environmental effects on mouse intestinal 46. Sun X, Shen J, Liu C, Li S, Peng Y, Chen C, et al. L-Arginine and N-carba‑ microbiota. ISME J. 2012;6:2033–44. moylglutamic acid supplementation enhance young rabbit growth and 21. Cartuche L, Pascual M, Gómez EA, Blasco A. Economic weights in rabbit immunity by regulating intestinal microbial community. Asian Australas J meat production. World Rabbit Sci. 2014;22:165–77. Anim Sci. 2019;33:166–76. 22. Zeng B, Han S, Wang P, Wen B, Jian W, Guo W, et al. The bacterial commu‑ 47. Piles M, Sánchez JP. Use of group records of feed intake to select for feed nities associated with fecal types and body weight of rex rabbits. Sci Rep. efficiency in rabbit. J Anim Breed Genet. 2019;136:474–83. 2015;5:9342. 48. de Agüero MG, Ganal-Vonarburg SC, Fuhrer T, Rupp S, Uchimura Y, Li 23. Drouilhet L, Achard CS, Zemb O, Molette C, Gidenne T, Larzul C, et al. H, et al. The maternal microbiota drives early postnatal innate immune Direct and correlated responses to selection in two lines of rabbits development. Science. 2016;351:1296–302. selected for feed efficiency under ad libitum and restricted feed‑ 49. Lagkouvardos I, Lesker TR, Hitch TC, Gálvez EJ, Smit N, Neuhaus K, et al. ing: I. Production traits and gut microbiota characteristics. J Anim Sci. Sequence and cultivation study of Muribaculaceae reveals novel species, 2016;94:38–48. host preference, and functional potential of this yet undescribed family. 24. Velasco-Galilea M, Piles M, Ramayo-Caldas Y, Sánchez JP. The value of gut Microbiome. 2019;7:28. microbiota to predict feed efficiency and growth of rabbits under differ‑ 50. Pereira FC, Wasmund K, Cobankovic I, Jehmlich N, Herbold CW, Lee KS, ent feeding regimes. Sci Rep. 2021;11:19495. et al. Rational design of a microbial consortium of mucosal sugar utilizers 25. Xu L, Paterson AD, Turpin W, Xu W. Assessment and selection of reduces Clostridiodes difficile colonization. Nat Commun. 2020;11:5104. competing models for zero-inflated microbiome data. PLoS One. 51. Goodrich JK, Davenport ER, Beaumont M, Jackson MA, Knight R, Ober C, 2015;10:e0129606. et al. Genetic determinants of the gut microbiome in UK twins. Cell Host 26. Gómez EA, Rafel O, Ramón J. The caldes strain (Spain). Options Mediterr B. Microbe. 2016;19:731–43. 2002;38:193–8. 52. Bonder MJ, Kurilshikov A, Tigchelaar EF, Mujagic Z, Imhann F, Vila AV, 27. Desjardins P, Conklin D. NanoDrop microvolume quantitation of nucleic et al. The effect of host genetics on the gut microbiome. Nat Genet. acids. J Vis Exp. 2010;45:2565. 2016;48:1407–12. 28. Parada AE, Needham DM, Fuhrman JA. Every base matters: assess‑ 53. Abbas W, Howard JT, Paz HA, Hales KE, Wells JE, Kuehn LA, et al. Influence ing small subunit rRNA primers for marine microbiomes with mock of host genetics in shaping the rumen bacterial community in beef cat‑ communities, time series and global field samples. Environ Microbiol. tle. Sci Rep. 2020;10:15101. 2016;18:1403–14. 54. Le Roy CI, Beaumont M, Jackson MA, Steves CJ, Spector TD, Bell JT. Herit‑ 29. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello able components of the human fecal microbiome are associated with EK, et al. QIIME allows analysis of high-throughput community sequenc‑ visceral fat. Gut Microbes. 2018;9:61–7. ing data. Nat Methods. 2010;7:335–6. 55. Xu F, Fu Y, Sun TY, Jiang Z, Miao Z, Shuai M, et al. The interplay between 30. Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME host genetics and the gut microbiome reveals common and distinct improves sensitivity and speed of chimera detection. Bioinformatics. microbiome features for complex human diseases. Microbiome. 2011;27:2194–200. 2020;8:145. 31. Rideout JR, He Y, Navas-Molina JA, Walters WA, Ursell LK, Gibbons SM, 56. La Reau AJ, Meier-Kolthoff JP, Suen G. Sequence-based analysis of the et al. Subsampled open-reference clustering creates consistent, com‑ genus Ruminococcus resolves its phylogeny and reveals strong host prehensive OTU definitions and scales to billions of sequences. PeerJ. association. Microb Genom. 2016;2:e000099. 2014;2:e545. 57. Gonzalez A, Navas-Molina JA, Kosciolek T, McDonald D, Vázquez-Baeza Y, 32. McDonald D, Price MN, Goodrich J, Nawrocki EP, DeSantis TZ, Probst Ackermann G, et al. Qiita: rapid, web-enabled microbiome meta-analysis. A, et al. An improved Greengenes taxonomy with explicit ranks for Nat Methods. 2018;15:796–8. ecological and evolutionary analyses of bacteria and archaea. ISME J. 2012;6:610–8. 33. Paulson JN, Stine OC, Bravo HC, Pop M. Differential abundance analysis Publisher’s Note for microbial marker-gene surveys. Nat Methods. 2013;10:1200–2. Springer Nature remains neutral with regard to jurisdictional claims in pub‑ 34. Henderson CR. Sire evaluation and genetic trends. J Anim Sci. lished maps and institutional affiliations. 1973;1973(Issue Symposium):10–41. 35. Rodrigues-Motta M, Gianola D, Heringstad B, Rosa GJM, Chang YM. A zero-inflated Poisson model for genetic analysis of the number of mastitis cases in Norwegian Red cows. J Dairy Sci. 2007;90:5306–15. 36. Varona L, García-Cortés LA, Pérez-Enciso M. Bayes factors for detection of quantitative trait loci. Genet Sel Evol. 2001;33:133–52.