Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Genetics Selection Evolution https://doi.org/10.1186/s12711-023-00803-4 RESEARCH ARTICLE Open Access On the holobiont ‘predictome’ of immunocompetence in pigs Joan Calle‑García1, Yuliaxis Ramayo‑Caldas2, Laura M. Zingaretti1, Raquel Quintanilla2, María Ballester2 and Miguel Pérez‑Enciso1,3,4*    Abstract Background Gut microbial composition plays an important role in numerous traits, including immune response. Integration of host genomic information with microbiome data is a natural step in the prediction of complex traits, although methods to optimize this are still largely unexplored. In this paper, we assess the impact of different model‑ ling strategies on the predictive capacity for six porcine immunocompetence traits when both genotype and micro‑ biota data are available. Methods We used phenotypic data on six immunity traits and the relative abundance of gut bacterial communities on 400 Duroc pigs that were genotyped for 70 k SNPs. We compared the predictive accuracy, defined as the correla‑ tion between predicted and observed phenotypes, of a wide catalogue of models: reproducing kernel Hilbert space (RKHS), Bayes C, and an ensemble method, using a range of priors and microbial clustering strategies. Combined (holobiont) models that include both genotype and microbiome data were compared with partial models that use one source of variation only. Results Overall, holobiont models performed better than partial models. Host genotype was especially relevant for predicting adaptive immunity traits (i.e., concentration of immunoglobulins M and G), whereas microbial composi‑ tion was important for predicting innate immunity traits (i.e., concentration of haptoglobin and C-reactive protein and lymphocyte phagocytic capacity). None of the models was uniformly best across all traits. We observed a greater variability in predictive accuracies across models when microbiability (the variance explained by the microbiome) was high. Clustering microbial abundances did not necessarily increase predictive accuracy. Conclusions Gut microbiota information is useful for predicting immunocompetence traits, especially those related to innate immunity. Modelling microbiome abundances deserves special attention when microbiability is high. Clus‑ tering microbial data for prediction is not recommended by default. Background Immunocompetence is an important contributor to pro- ductivity, profitability, and welfare in pigs as in other spe- *Correspondence: Miguel Pérez‑Enciso cies [1, 2]. The current emergence of antibiotic resistance mperezenciso@gmail.com [3], coupled with increasing social demands for healthier 1 Centre for Research in Agricultural Genomics CSIC-IRTA-UAB-UB, products and environmentally responsible livestock sys- Campus UAB, Edifici CRAG​, 08193 Bellaterra, Spain 2 Animal Breeding and Genetics Program, Institut de Recerca i Tecnologia tems, favour the incorporation of health-related traits in Agroalimentàries (IRTA), Caldes de Montbui, 08140 Barcelona, Spain pig breeding programs [4]. Importantly, the relevance of 3 ICREA, Passeig Lluis Companys 23, 08010 Barcelona, Spain the pig as a biomedical model is rapidly increasing, given 4 Present Address: Corteva Agriscience, Virtual Location, Bergen op Zoom, Indianapolis 4611 BB, Netherlands its similarities with humans in terms of anatomy and physiology [5, 6]. © The Author(s) 2023. 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. Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 2 of 12 The importance of the composition of human and ani- immunocompetence in pigs. The diversity of the genetic mal microbial symbionts for health and production is and microbial influences on immunocompetence in pigs now widely recognized [7–11]. In this context, the ‘hol- that has been reported in our previous work [9] makes ogenome’ concept [12], which describes the joint action immunocompetence traits a good case study for explor- of host genome and microbiome on a phenotype, is ing the impact of different holobiont modelling strate- becoming increasingly popular in the literature. Since gies. Here, we focus on prediction and explore a wide the pioneering work of Ross et al. [13], in which human range of choices in terms of statistical models, priors, and health-related traits and methane production of a cow abundance clustering. We named this large catalogue of were predicted using microbiota information, numer- model choices for prediction the ‘predictome’. ous studies have used the composition of gut microbial ecosystems to elucidate its contribution to different com- Methods plex animal traits, such as methane emission in cattle Animal samples [13–15], carcass traits in pigs [16, 17], and feed efficiency The animal material used in this study belonged to a in several farm animals [18–22], among others. The rel- commercial Duroc line and is fully described in Ramayo- evance of gut microbiota composition in the mammalian Caldas et al. [9]. The animals used here were a subset of immune system has received special attention in the past 400 Duroc piglets (199 females and 201 males) out of the years [10, 23]. Several studies have reported associations 432 in [9]. They were the offspring of 22 boars and 132 between microbial abundances and genomic polymor- sows. Animals were raised on the same farm, fed ad libi- phisms in immunity-related genes [24–27], which is in tum with a commercial cereal-based diet. They were line with previous reports stating that the host genome apparently healthy, without any sign of infection. Blood can influence the composition of the gut microbiota and faecal samples were collected at 60 ± 8 days of age in [26, 28–30]. In a recent study [9], we explored the con- six different batches (dates). Fecal samples were trans- tribution of the host’s genotype and its gut microbiota to ferred to cryotubes, conserved in ice, and later stored several immunocompetence traits in a Duroc pig popula- at − 80 °C until DNA extraction. tion, and we found an important connection between gut microbiota composition and pig immunity. Immunocompetence traits Considering the contribution of microbial composition For this study, we selected six phenotypes that covered to animal physiology and immunity, the next natural step a range of genetic and microbial parameters based on is to incorporate this information into a prediction frame- previous evidence [9]. These phenotypes were classified work. However, guidelines for the integration of micro- according to the two major components of the immune bial abundances with the currently ubiquitous genotype system: plasma concentrations of immunoglobulins M marker data are lacking. The analogy between microbi- (IgM) and G (IgG), which are associated with adaptive ome and genotype marker data, i.e., that the number of immunity, and serum concentrations of the acute-phase features typically exceeds the number of samples, should proteins haptoglobin (HP) and C-reactive protein (CRP) not hide some important differences between both types plus the phagocytic capacity of lymphocytes (LYM_ of data: (i) microbiota composition may change through- PHAGO_FITC), which are components of innate immu- out the life of an organism [31], (ii) there are many more nity. Finally, the percentage of gamma-delta T cells (γδ microbial genes than genes in the host genome [32], (iii) T cell), which are a bridge between innate and adaptive the microbiome has a hierarchical structure given by its immunity, was also considered. Details of the sampling phylogenetic relationships, (iv) microbial abundance and laboratory processing are in [9]. All traits, except distributions are zero-inflated and highly leptokurtic, immunoglobulins (IgM and IgG) had a highly leptokur- and (v) microbial abundances are compositional, which tic distribution, which was addressed by log-transfor- induces relationships between abundances [33]. mation. Data were preadjusted prior to cross-validation Using the holobiont, i.e., genotype and microbiome to focus on genetic and microbial effects rather than on data, to better understand and predict complex traits is environmental factors. Following previous analyses [39], still largely unexplored. Some studies have integrated IgM, IgG, CRP, and γδ T cell phenotypes were corrected genotype and microbiota composition data in the same for the effect of batch (sampling date), HP was corrected model to explore their contributions to complex traits for the effects of batch and sex, and LYM_PHAGO_FITC [18, 34–36], while genetic and microbial parameters were was corrected for the effects of batch and sex. independently analysed in other studies [15, 21, 22, 37]. In this study, and complementing our previous sim- Genotype data ulation work [38], we evaluated the performance of Animals were genotyped for 68,516 single nucleotide pol- alternative modelling strategies for the prediction of ymorphisms (SNPs) with the Porcine 70 k GGP Porcine Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 3 of 12 HD Array (Illumina, San Diego, CA, USA), using the the number of variables in the model. Two clustering Infinium HD Assay Ultra protocol (Illumina). We used options were tested: by phylogeny and by abundance. the Plink v1.9 software [40, 41] to remove SNPs that had For taxonomic assignment, a classifier was created a minor allele frequency less than 5%, that had more than with QIIME2 v 2021.8 using the GreenGenes version 10% missing genotype data, that mapped to the sex chro- 13.8 database [48]. For phylogeny clustering, the ASV mosome, or that did not map to the porcine reference sequences that passed quality control were aligned using genome (Sscrofa11.1 assembly). After quality control, the MUSCLE algorithm [49] and the output was used to 41,131 SNPs were retained for subsequent analysis. Miss- perform the phylogeny analysis with the UPGMA algo- ing genotypes were rare (0.19%) and were simply imputed rithm [50]. Alignment and phylogeny analyses were con- with the average allele frequency of each SNP. ducted using MEGA version 11 [51]. The resulting rooted phylogeny tree was processed with the cutree function Microbiota abundance data of the stats R package [52]. Since not all ASV could be To increase the read depth of microbial abundances assigned to the genus level and not all genuses collapsed available in Ramayo-Caldas et al. [9] and increase the at the same height of the phylogeny tree, the tree was reliability of the experiment reported here, we combined sliced at different heights (h), and the consistency in the 16S metagenome data from [9] with new sequence data taxonomic composition at the genus level of the ASV in obtained from the same DNA samples. The bioinformatic each subtree was evaluated. We found height h = 0.6 to procedure was the same and a strict quality control was be the optimum value. The leaves, i.e., the ASV, of each applied to ensure the data were comparable, as described subtree were clustered, yielding k = 232 phylogeny-based below. clusters that corresponded approximately to the genus DNA extraction was carried out with the DNeasy Pow- levels. The hclust package with Ward’s [53] method was erSoil Kit (QIAGEN, Hilden, Germany), following manu- used for abundance clustering. The cutree function was facturer’s instructions. The first sequence dataset was used again to generate the number of desired clusters. obtained with paired end NovaSeq (2 × 250 nucleotides), We used k = 232 to obtain the same number of clusters while the second sequence batch was obtained with as in the phylogeny approach and facilitate comparison paired end Illumina MiSeq (2 × 300). The same primers between clustering methods. were used for both batches: (V3_F357_N) 5′-CCT​ACG​ GGNGGC​WGC​AG-3′ and (V4_R805) 5′-GAC​TAC​ Prediction analysis HVGGG​TAT​CTA​ATC​C-3′. We compared the predictive accuracy of Bayesian repro- Each 16S sequence dataset was processed indepen- ducing kernel Hilbert space (RKHS) regression and Bayes dently with the QIIME2 v 2021.8 software [42], using C [54] models implemented with the BGLR R pack- the same bioinformatic pipelines. Denoising [43–46] was age [55]. In Bayes C, the combined, holobiont model performed to extract each amplicon sequence variant (referred to as model XB) was: (ASV) from the raw sequencing data using the R pack- age dada2 [44], as implemented in the denoise-paired y = m + XβG + BβB + e, (1a) QIIME2 plugin. Primers were manually removed in both where y is the corrected phenotype vector of size n, batches by trimming out the first 17 and 21 nucleotides the number of samples (400); m is the general mean; X from the forward and reverse strands, respectively. No is the matrix of standardized SNP genotype value; B is truncation was performed on the first sequencing batch, the standardized CLR-transformed matrix of microbial but both forward and reverse sequences were truncated abundances, either individual ASV or clusters; βG and βB to 250 nucleotides in the second batch due to low qual- are the corresponding vectors of SNP and ASV effects, ity. The two batches were merged after denoising with the respectively; and e is the vector of residuals. Prior to feature-table merge plugin and the overlap method ‘sum’. standardization, SNPs were originally coded as 0, 1, 2 for The raw number of ASV was 57,195, but ASV that were ‘AA’, ‘AB’ and ‘BB’ genotypes, respectively. Partial models present in less than three samples and that represented that considered only genotypes (model X) or ASV abun- less than 0.001% of the total counts were discarded. These dances (model B) were also evaluated: cut-offs are within the optimal range suggested by [47]. Centered log-ratio (CLR) transformation was applied to y = m + XβG + e, (1b) raw ASV abundances for further analyses. The pipeline used to process and combine these datasets is in Addi- y = m + BβB + e. (1c) tional file 1. Abundances were considered either by ASV or clus- In the case of the kernel-based method RKHS, the hol- tered, the latter in order to evaluate the effect of reducing obiont and partial models were: Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 4 of 12 Estimation of heritability and microbiability y = m + Zu + Wb + e, (2a) We estimated heritability ­ (h2) and microbiability ­ (b2) using all combinations of priors and models described, y = m + Zu + e, (2b) using now the complete dataset. For RKHS, variance components are explicitly defined. For the Bayes C mod- y = m + Wb + e, (2c) els, we used the approach suggested in [55] (https://​ github.​ com/​ g dlc/​ B GLR-R/​ b lob/​ m aster/​ i nst/​ m d/​ h erit​ respectively, where Z and W are incidence (here iden- abili​ty.​md) to estimate heritability and microbiability. In tity) matrices for genotype and microbial abundances, short, at each iteration i, the method samples the effects respectively; u and b are vectors with genotypic and of SNPs and ASV: microbial random effects; and e is the vector of residuals. We assumed u ∼ N (0, Gσg2 ) and b ∼ N (0, Mσb2 ), where u(i) = XβG (i) , (3) XX ′ BB′ G= nG and M = nBare the genomic and microbial relationship matrices, respectively, nG being the number b(i) = BβB(i) , (4) of SNPs (41,131) and nB, the number of ASV (2945). We compared a range of priors for both abundances where u(i) and b(i) are sampled genome and micro- and genotypes. In Bayes C, we assessed varying prior biota effects at the i -th iteration for the set of individu- probabilities π0 of a feature (SNP, ASV or cluster) to als, respectively. Therefore, h(i)2 = Var(u )/Var(y) and (i) enter the model. For SNPs, the values used were π0 = 0.5, b(i) = Var(b(i) )/Var(y) are the sampled h 2 ­ 2 and ­b2 in the 0,01, 0.001, and 0.0001; for individual ASV, π0 = 0.5, 0.1, i -th iterate, from which posterior means were estimated 0.01, and 0.001; and for ASV clusters, π0 = 0.5, 0.1 and by averaging over iterations. 0.01. In the case of RKHS, default prior parameters for In addition, we computed the contribution of each ASV variance components were taken as informative (referred to total microbiability by setting the estimated effects for to as FALSE), whereas an uninformative REML-like prior all ASV to zero except the ASV of interest. This was done (referred to as TRUE) was achieved with parameters using the modelling combination that yielded the high- df0 = 0.0001 and S0 = 0.0001, where df0 is the degrees est predictive accuracy and that used Bayes C to model of freedom and S0 is the scale parameter of the variance microbial abundances. We also estimated the heritability component. of the abundance of each ASV using RKHS. BGLR allows each feature, genotypes or abundances, to be modelled independently; thus we considered several modelling combinations: the mentioned methods and Results Quality control priors under both partial and combined models. Abun- Several quality control measures were applied to ensure dances were considered either at the individual ASV level that 16S reads from the two sequencing batches could or clustered. In total, 133 different models were evaluated be merged: (i) Euclidean distances between samples for for each trait. In addition, we also compared an ensemble each pair of datasets were highly correlated (Additional method, computed as the average predicted values from file 3: Table S1), (ii) samples showed no structure in a all analyses for that trait [56]. principal component analysis (PCA), neither for individ- To evaluate the accuracy of predictions, three parti- ual datasets nor for the merged dataset (see Additional tions of 80 (20%) samples were removed and predicted file 2: Fig. S2) and (iii) at least 76% of ASV detected in using the model that was trained with the remaining the post-denoising combined dataset were also detected data. Partitions were randomly chosen but keeping a in each individual batch, and the number of detected fixed proportion of samples within each of the sampling ASV was similar for the two batches: 2971 vs. 2566 for batches. The correlation between observed and predicted batches 1 and 2, respectively. In the merged dataset, 2945 phenotypes, averaged over the three partitions, was used ASV pertaining to 53 genera were detected, of which 54% as a measure of predictive accuracy. In each analysis, are classified at the genus level. In this dataset, the aver- BGLR was run for 100k iterations, including 500 burn-in age number of 16S reads per sample was 136,616 (Addi- iterations, and thinning every five iterations. This num- tional file 2: Fig. S3). As expected, the distribution of ber of iterations seemed sufficient for convergence (see ASV frequencies was highly leptokurtic-over 75% of the Additional file 2: Fig. S1). Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 5 of 12 Fig. 1 Predictive accuracy for each model and prior combination considered. Predictive accuracy, defined as correlation between predicted and observed phenotypes for each model, averaged over the three partitions. Models are classified by data used: SNPs only (Model X, Eqs. (1b) and (2b)), ASV only (Model B, Eqs. (1c) and (2c)), both SNPs and ASV (Model XB, Eqs. (1a) and (2a)). Each dot corresponds to the prediction accuracy obtained with a different combination of statistical method and prior, averaged over the three partitions microbial ASV were present in only 124 samples or less all traits. However, it is relevant to note that the ensem- (Additional file 2: Fig. S4). ble method was consistently better, or very similar in the worst case, than the average performance of the two indi- Predictive accuracies vidual methods. Beyond that, no clear pattern emerged. Predictive accuracies across models are shown in Fig. 1. Given that abundance distributions were highly lep- Each dot corresponds to a predictive accuracy obtained tokurtic (Additional file 2: Fig. S4), it is pertinent to ask with a different combination of statistical method and whether clustering ASV in fewer groups could improve prior, averaged over the three partitions (see “Methods”). prediction. Figure 3 shows prediction accuracies when The relevance of genotypes and microbiota in predic- all 2945 ASV were included in the model individually, tion varied. We found genotype information to be more clustered at the genus level or by abundance. For the relevant than microbial abundances for predicting IgM latter two cases, the number of clusters was k = 232. and IgG and for the proportion of γδ T cells. In con- Only traits for which microbial abundances were rel- trast, genotype data was not relevant for predicting HP. evant for prediction (Fig. 1) were considered in this For CRP, both sources of information improved predic- comparison, i.e., IgM, HP, CRP, LYM_PHAGO_FITC. tion when considered jointly: the best partial models X Clustering had an important effect on prediction accu- (Eqs. (1b) and (2b)) and B (Eqs. (1c) and (2c)) yielded racy for all traits, although not always in a positive similar predictive accuracies, which notably increased direction. Clustering by abundance improved predic- with the combined XB models (Eqs. (1a) and (2a)). Nev- tion for HP and by phylogeny for CRP concentration. ertheless, the best combined model was the best strategy Clustering worsened prediction for the remaining for all traits, except for a single analysis under a partial traits, especially for LYM_PHAGO_FITC. microbiota model (B) for LYM_PHAGO_FITC. However, far more interesting is the observation that variability in predictions, i.e., ‘sensitivity’ to modelling, of models that Impact of modelling on heritability and microbiability include microbiota information was larger for traits for estimates which the microbiome was more relevant (HP, CRP and Different modelling approaches are expected to be LYM_PHAGO_FITC) than for traits for which genotypes reflected in estimates of ­h2 and b ­ 2. Furthermore, although sufficed for prediction (immunoglobulins and γδ T cells). the ‘true’ model is never known, it is useful to know how Next, we investigated the impact of RKHS and Bayes C sensitive estimates of ­h2 and ­b2 are to alternative mod- on prediction using different combinations of holobiont els, and what are the estimates that correspond to the models (Fig. 2). No method was uniformly best across best models in terms of prediction. Unsurprisingly, given Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 6 of 12 Fig. 2 Predictive accuracy for each model in holobiont models. Predictive accuracy, defined as correlation between predicted and observed phenotypes for each model, averaged over the three partitions. The first element of each combination is the method (Bayes C or RKHS) used to model the genotype, while the element after the dot is the method used to model microbiota abundance. Each dot corresponds to the prediction accuracy obtained with a different combination of statistical method and prior, averaged over the three partitions. The dashed line is the predictive accuracy of the ensemble method Fig. 3 Predictive accuracy with microbiome clustering for the partial microbiota model. Predictive accuracy, defined as correlation between predicted and observed phenotypes for each analysis and trait, averaged over the three partitions, for each clustering approach: NONE, microbiota at the ASV level; ABUND, microbiota clustered by abundance; PHYLO, microbiota clustered by phylogeny. Results for IgG and γδ T cells are not shown since their prediction accuracy was not influenced by the microbiome (Fig. 1). Each dot corresponds to the prediction accuracy obtained with a different combination of statistical method and prior, averaged over the three partitions the prediction results (Fig. 1), estimates of ­h2 were much especially for ­b2. Overall, larger estimates of either ­h2 larger than estimates of ­b2 for IgG and γδ T cells, while or ­b2 were obtained with the simple X or B models than the opposite was observed for HP (Fig. 4). Estimates of with the complete holobiont models. This trend was far ­h2 and b­ 2 obtained from the holobiont (XB) model were more marked for the estimates of b ­ 2 and suggests some comparable in the case of CRP. We observed a large confounding between genotype and microbiome effects. variability in estimates across modelling combinations, Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 7 of 12 Fig. 4 Distribution of estimates of heritability (a) and microbiability (b) obtained from partial and combined models. Estimates were obtained from the whole dataset. Model X: only genotypes are used (Eqs. (1b) and (2b)); model B: only abundances are used (Eqs. (1c) and (2c)); model XB: genotypes and abundances used (Eqs. (1a) and (2a)). Analyses with clustered microbial abundances were excluded. Each dot corresponds to an estimate obtained with a different combination of statistical method and prior in the whole dataset Effect of microbial abundances on immunocompetence the phenotypic variation in LYM_PHAGO_FITC. It is Figure 5 shows the cumulative sum of the contributions relevant to note that some of these ASV exhibited much of ASV to the phenotypic variance (b2j), inferred from a higher ­h2 estimates than average (see Additional file 2: Bayes C analysis. Curves for IgG and percentage of γδ T Fig. S5). Several ASV that were annotated at the genera cells clearly indicated that the contribution of the micro- level had particularly high ­h2 estimates (Treponema sp, biome was negligible for these traits. As suggested by ­h2 ~ 0.6; Streptococcus, ­h2 ~ 0.5, and (see Additional file 3: the straight line, estimates of the effects of ASV on the Table S2). phenotypes were uniformly distributed for CRP and IgM. For HP, a few ASV had somewhat larger effects, but esti- Discussion mated effects were similar across the other ASV. At the Is it worth using holobiont data to improve prediction? opposite extreme, the cumulative distribution of esti- In a previous simulation study [38], we hypothesized mates of effects on LYM_PHAGO_FITC indicated that that microbiome data might increase predictive ability a few ASV explained a large part of b ­ 2. For this trait, 19 of complex phenotypes by up to ~ 50% in an optimistic 2 ASV were responsible for 30% of b ­ (see Additional file 3: scenario where ­h2 ~ ­b2 ~ 0.25 with stable microbiota com- Table S2), i.e., 10% of the total phenotypic variation of position and large sample size (~ 1k). In general, the lit- this trait. erature reports a positive impact of combining genotype The effects of the few ASV that were responsible for a and microbiome data for the prediction of phenotype large proportion of ­b2 for LYM_PHAGO_FITC can rep- for complex traits, although results vary according to resent direct or indirect genetic effects. Given the special trait and species (Table 1). The results presented here are microbial architecture of this trait (Fig. 5), we explored broadly consistent with these studies. the heritability of the 19 ASV with the largest contribu- The advantage of adding microbiota for prediction, tion (see Additional file 3: Table S2). ASV belonging to on top of genotype data, will depend on its stability and the Bacteroidales order and to the Ruminococcaceae on the ratio of microbiability to heritability. Ramayo- family together explained ~ 5% of the phenotypic vari- Caldas et al. [9] reported higher estimates of h ­ 2 than of ance in LYM_PHAGO_FITC; Treponema genera and ­b2 for concentrations of immunoglobulins G and M and Lactobacillales family were jointly responsible for 3% of the proportion of γδ T cells in pigs of the same age, while Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 8 of 12 Fig. 5 Cumulative sum of the contribution of individual ASV to microbiability for each trait. ASV microbiabilities from the modelling combination that yielded the highest prediction accuracy and that used Bayes C to model microbial abundances Table 1 Summary of the results in the literature for the prediction of the phenotype for complex traits using genotype and/or microbiome data Species Phenotypes N Host genotype Microbiome Results Refs. Cow Methane emissions; rumen and blood metabolites; milk 1016 120k SNPs 512 OTU B>X [15] production efficiency Cow Milk fatty acid content 292 39k SNPs 3055 OTU XB ~ X [35] Cow Milk acetone and β-hydroxybutyric acid 277 436k SNPs 4226 OTU B>X [37] Rabbit Feed efficiency; growth 425 Pedigree 963 OTU XB > X [18] Japanese quail Feed intake; P utilization; body weight gain; feed efficiency 758 4k SNPs 1188 OTU B>X [22] Pig Daily gain; feed conversion; feed intake 207 52k SNPs 2714 OTU B>X [21] Pig Meat quality 1123 42k SNPs 1755 OTU XB > X [34] Human Dietary habits; anthropometric and blood measurements 715 545k SNPs 1.36 M genes XB > X [36] (metagenomics) N: Sample size B > X: performance using microbiome was larger than with genotypes only; XB > X: holobiont model was better than using genotypes only; XB ~ X: both models performed similarly estimates of ­h2 and ­b2 were similar for traits related to outperformed, or as least matched, the best partial model, innate immunity (HP, CRP, and phagocytic capacity of except for LYM_PHAGO_FITC (Fig. 1). In the case of HP lymphocytes). Here we show that, by default, the safest concentration, for which the estimate of microbiability choice is the use of a holobiont model rather than partial was high, the predictive accuracy increased when micro- X or B models. We observe that the best holobiont model biome information was added (0.29), compared to using Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 9 of 12 only genotype data (0.04). For other traits, such as IgM correlation may arise because they are evolutionarily and CRP, the increase in predictive accuracy for the XB related and can therefore share the same ecological niche vs the X model was ~ 10 and 30%, respectively. Therefore, [59], or can form ‘micro communities’ for which abun- although microbiota data generally improves prediction dances of several microorganisms follow similar pat- and a holobiont model is recommended, gains can be terns [60]. Thus, we hypothesized that clustering ASV limited compared to the best-case simulated scenario in abundances could result in improved predictions, given [38], which was up to 50% that the abundance of clusters would have better distri- butional properties and would average out the noise of Impact of the model used for prediction the abundance of individual ASV. Clustering did affect While using DNA markers for complex trait prediction predictions but not always positively. We found that clus- has a long history and several statistical models have tering improved prediction in cases when the best micro- become widely used (e.g., genomic best linear unbiased biome modelling option was RKHS, such as for HP and prediction (GBLUP)), much less is known on optimum CRP. In contrast, clustering worsened predictions for methods to leverage microbiome data. A major reason LYM_PHAGO_FITC, where Bayes C was the best model. is that distributions of microbiome abundance are highly Our hypothesis is that clustering can mask strong micro- leptokurtic and remain non-normal even after usual bial signals when there are few causative unrelated bacte- transformations [33]. In addition, stability of microbial ria by grouping them with irrelevant ASV. Nevertheless, composition is controversial, especially for rare microbes taking each ASV individually seems the safest default [31]. We and others have argued that there is much room choice. for methodological improvement in this area [38, 63]. We observed that ‘sensitivity’ to modelling, i.e., the Dissecting the influence of microbiota on complex traits variability in predictive accuracies across modelling com- In our previous simulation work [38], we showed that binations, differed markedly between traits. Interestingly, distinguishing between direct and indirect genetic this sensitivity increased with microbiability. This can effects mediated by the microbiota on a complex phe- be appreciated in Fig. 1 by comparing the dispersion of notype is difficult. Here, estimates of the ­h2 were low combined model predictive accuracies for traits related for the abundance of most ASV (see Additional file 2: to adaptative immunity (IgM, IgG and γδ T cells) vs. the Fig. S5, median = 0.07). This is consistent with previous innate immunity traits (HP, CRP and LYM_PHAGO_ studies on the genetic control of host genotype over gut FITC). Then, unsurprisingly, accuracies of models that microbial composition in pigs [21, 61] and humans [36], included abundances were more variable than those that and suggests that indirect genetic ASV-mediated effects included only genotypes. All this indicates that modelling should be rare or may not have a large impact on complex of microbial abundances deserves special care. This was phenotypes. However, and as previously reported [18, observed irrespective of whether Bayes C or RKHS was 34–36], there are traits for which putative causative ASV used, which suggests an important role of prior infor- can in turn be highly heritable. For LYM_PHAGO_FITC, mation. Interestingly, for all traits, ensemble method the abundance of a few ASV explained a large part of the predictions were consistently better, or at worse very variability in that trait (Fig. 5), and some of them exhib- similar, than the average prediction obtained for almost ited high estimates of h ­ 2 (see Additional file 3: Table S2). all modelling combinations. Therefore, ensemble meth- In this case, a partial indirect model cannot be ruled out. ods provide a safe and ‘agnostic’, albeit computationally Bacteria members of the Spirochaetaceae, Prevotellaceae expensive, approach for prediction—and at the cost of and Streptococcaceae families were the main candidates hindering biological interpretation [57]. for explaining putative indirect effects, given their con- No model was uniformly best across all traits. In our tribution to the trait variability and their high heritability previous simulation study [38], we observed that Bayes C estimates. Previous studies also reported high estimates tended to predict phenotype better than RKHS, although of ­h2 for Prevotellaceae and Streptococcaceae in the the advantage decreased as the number of causative ASV swine gut [61], and for Prevotellaceae in the human gut increased. Here, we found that Bayes C was the best [62]. However, this speculation must be considered with model for three traits: IgG, CRP, and LYM_PHAGO_ care, given the low prevalence of some ASV and the high FITC. This suggests that the number of relevant features error level of individual estimates of heritability. Taken is smaller for these three traits than for the other traits, together, our results indicate that more sophisticated sta- for which an infinitesimal model would fit the data better tistical analyses are needed to get a better understanding [58]. of how direct and microbiome-mediate genetic effects Microbial abundances do not only have a skewed dis- can be dissected [11, 63]. tribution, they are also correlated between them. This Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 10 of 12 Conclusions Availability of data and materials The raw sequencing data of host genotype and both batches of 16S gut Microbiota information improves the prediction of metagenome employed in this article are available in the NCBI’s short read phenotypes for immunocompetence traits, especially archive (https://​www.​ncbi.​nlm.​nih.​gov/​sra) BioProject PRJNA608629. for traits related to innate immunity. As a result, the best default option is a holobiont model, rather than Declarations partial models that use only genotype or microbial data. Ethics approval and consent to participate Care in choosing the statistical model is particularly Animal care and experimental procedures were carried out following national important for traits for which microbiability is high. and institutional guidelines for the Good Experimental Practices and were Somewhat counter-intuitively, clustering microbial taxa approved by the IRTA Ethical Committee. Consent to participate is not applica‑ ble in this study. does not necessarily help prediction, but if it does, it is for traits for which most microbial ASV have similar, Consent for publication small effects. The best clustering strategy is trait-spe- Not applicable. cific. Confounding of the estimates of heritability and Competing interests microbiability for a trait may suggest an indirect effect The authors declare that they have no competing interests. of host genotype through microbial composition. Evi- dence in favour of indirect genetic effects seems weak Received: 5 October 2022 Accepted: 7 April 2023 though, more sophisticated statistical models than those used here are warranted to settle this issue. Supplementary Information References The online version contains supplementary material available at https://​doi.​ 1. Heuß EM, Pröll-Cornelissen MJ, Neuhoff C, Tholen E, Große-Brinkhaus org/​10.​1186/​s12711-​023-​00803-4. C. Invited review: piglet survival: benefits of the immunocompetence. Animal. 2019;13:2114–24. 2. Aleri JW, Hine BC, Pyman MF, Mansell PD, Wales WJ, Mallard B, et al. Additional file 1: QIIME2 pipeline used to process and combine the two Associations between immune competence, stress responsiveness, and 16S sequence datasets. production in Holstein-Friesian and Holstein-Friesian × Jersey heifers Additional file 2: Fig. S1. Evolution of heritability estimates along itera‑ reared in a pasture-based production system in Australia. J Dairy Sci. tions in an holobiont model in CRP to illustrate convergence of the MCMC 2019;102:3282–94. chain. Fig. S2. PCA for the individual batches and the merged dataset 3. Cheng G, Hao H, Xie S, Wang X, Dai M, Huang L, et al. Antibiotic alterna‑ of CLR-transformed ASV abundances. Fig. S3. Distribution of reads per tives: the substitution of antibiotics in animal husbandry? Front Microbiol. sample in the merged 16S sequencing data after quality control. Fig. 2014;5:217. S4. Distribution of CLR-transformed ASV abundances, averaged across 4. Reverter A, Hine BC, Porto-Neto L, Li Y, Duff CJ, Dominik S, et al. Immune‑ samples. Fig. S5. Distribution of the ASV heritability estimates. DEX: a strategy for the genetic improvement of immune competence in Additional file 3: Table S1. Correlation of the Euclidean distances Australian Angus cattle. J Anim Sci. 2021;99:skaa384. between samples in each pair of datasets. Table S2. ASV responsi‑ 5. Meurens F, Summerfield A, Nauwynck H, Saif L, Gerdts V. The pig: a model ble for 30% of the estimated b2 under the best predictive model for for human infectious diseases. Trends Microbiol. 2012;20:50–7. LYM_PHAGO_FITC. 6. Walters EM, Wells KD, Bryda EC, Schommer S, Prather RS. Swine models, genomic tools and services to enhance our understanding of human health and diseases. Lab Anim (NY). 2017;46:167–72. Acknowledgements 7. Ruff WE, Greiling TM, Kriegel MA. Host–microbiota interactions in The authors warmly thank all the technical staff from Selección Batallé S.A. for immune-mediated diseases. Nat Rev Microbiol. 2020;18:521–38. providing the animal material and their collaboration during the sampling. 8. Zhang Q, Difford G, Sahana G, Løvendahl P, Lassen J, Lund MS, et al. We also thank Jack Dekkers for the thorough and rigorous editing of the Bayesian modeling reveals host genetics associated with rumen manuscript. microbiota jointly influence methane emission in dairy cows. ISME J. 2020;14:2019–33. Author contributions 9. Ramayo-Caldas Y, Zingaretti LM, Pérez-Pascual D, Alexandre PA, Reverter MPE, YRC and LMZ conceived research. MB, YRC and RQ provided data. JCG, A, Dalmau A, et al. Leveraging host-genetics and gut microbiota to deter‑ MPE, YRC and LMZ performed research. MPE and JCG wrote the manuscript mine immunocompetence in pigs. Anim Microbiome. 2021;3:74. with help from the other authors. All authors read and approved the final 10. Lo BC, Chen GY, Núñez G, Caruso R. Gut microbiota and systemic immu‑ manuscript. nity in health and disease. Int Immunol. 2021;33:197–209. 11. Weishaar R, Wellmann R, Camarinha-Silva A, Rodehutscord M, Bennewitz Funding J. Selecting the hologenome to breed for an improved feed efficiency in Project funded by Ministry of Science and Innovation-State Research pigs—a novel selection index. J Anim Breed Genet. 2020;137:14–22. Agency (AEI, Spain, https://​doi.​org/​10.​13039/​50110​00110​33) grant number 12. Zilber-Rosenberg I, Rosenberg E. Role of microorganisms in the evolution PID2019-108829RB-I00 to MPE, AGL2016-75432-R and PID2020-112677RB-C21 of animals and plants: the hologenome theory of evolution. FEMS Micro‑ to MB and RQ, and AGL2017-88849-R to YRC. This work was also supported by biol Rev. 2008;32:723–35. grant CEX2019-000902-S funded by MCIN/AEI/https://​doi.​org/​10.​13039/​50110​ 13. Ross EM, Moate PJ, Marett LC, Cocks BG, Hayes BJ. Metagenomic predic‑ 00110​33 and by the CERCA Programme/Generalitat de Catalunya (Spain). tions: from microbiome to complex health and environmental pheno‑ JCG is funded by a PhD fellowship PRE2020-094436 from AEI. YRC is recipient types in humans and cattle. PLoS One. 2013;8:e73056. of a Ramón y Cajal post-doctoral fellowship (RYC2019-027244-I) funded by 14. Difford GF, Plichta DR, Løvendahl P, Lassen J, Noel SJ, Højberg O, et al. the Spanish Ministry of Science and Innovation. MB, RQ and YRC belong to a Host genetics and the rumen microbiome jointly associate with methane Consolidated Research Group AGAUR, ref. 2017SGR-1719. emissions in dairy cows. PLoS Genet. 2018;14: e1007580. Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 11 of 12 15. Wallace RJ, Sasson G, Garnsworthy PC, Tapio I, Gregson E, Bani P, et al. 37. Gebreyesus G, Difford GF, Buitenhuis B, Lassen J, Noel SJ, Højberg O, A heritable subset of the core rumen microbiome dictates dairy cow et al. Predictive ability of host genetics and rumen microbiome for productivity and emissions. Sci Adv. 2019;5:eaav8391. subclinical ketosis. J Dairy Sci. 2020;103:4557–69. 16. Khanal P, Maltecca C, Schwab C, Fix J, Tiezzi F. Microbiability of meat 38. Pérez-Enciso M, Zingaretti LM, Ramayo-Caldas Y, de los Campos G. quality and carcass composition traits in swine. J Anim Breed Genet. Opportunities and limits of combining microbiome and genome data 2021;138:223–36. for complex trait prediction. Genet Sel Evol. 2021;53:65. 17. Maltecca C, Lu D, Schillebeeckx C, McNulty NP, Schwab C, Shull C, et al. 39. Ballester M, Ramayo-Caldas Y, González-Rodríguez O, Pascual M, Predicting growth and carcass traits in swine using microbiome data Reixach J, Díaz M, et al. Genetic parameters and associated genomic and machine learning algorithms. Sci Rep. 2019;9:6574. regions for global immunocompetence and other health-related traits 18. Velasco-Galilea M, Piles M, Ramayo-Caldas Y, Sánchez JP. The value of in pigs. Sci Rep. 2020;10:18462. gut microbiota to predict feed efficiency and growth of rabbits under 40. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, different feeding regimes. Sci Rep. 2021;11:19495. et al. PLINK: A tool set for whole-genome association and population- 19. Wen C, Yan W, Mai C, Duan Z, Zheng J, Sun C, et al. Joint contributions based linkage analyses. Am J Hum Genet. 2007;81:559–75. of the gut microbiota and host genetics to feed efficiency in chickens. 41. Purcell S Chang C. PLINK v1.90b6.24. 2021. https://​www.​cog-​genom​ics.​ Microbiome. 2021;9:126. org/​plink/1.​9/ Accessed 27 Feb 2023. 20. Delgado B, Bach A, Guasch I, González C, Elcoso G, Pryce JE, et al. 42. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, Whole rumen metagenome sequencing allows classifying and predict‑ et al. Reproducible, interactive, scalable and extensible microbiome ing feed efficiency and intake levels in cattle. Sci Rep. 2019;9:11. data science using QIIME 2. Nat Biotechnol. 2019;37:852–7. 21. Camarinha-Silva A, Maushammer M, Wellmann R, Vital M, Preuss S, 43. Quince C, Lanzén A, Curtis TP, Davenport RJ, Hall N, Head IM, et al. Bennewitz J. Host genome influence on gut microbial composi‑ Accurate determination of microbial diversity from 454 pyrosequenc‑ tion and microbial prediction of complex traits in pigs. Genetics. ing data. Nat Methods. 2009;6:639–41. 2017;206:1637–44. 44. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. 22. Haas V, Vollmar S, Preuß S, Rodehutscord M, Camarinha-Silva A, Ben‑ DADA2: high-resolution sample inference from illumina amplicon data. newitz J. Composition of the ileum microbiota is a mediator between Nat Methods. 2016;13:581–3. the host genome and phosphorus utilization and other efficiency traits 45. Amir A, McDonald D, Navas-Molina JA, Kopylova E, Morton JT, Xu ZZ, in Japanese quail (Coturnix japonica). Genet Sel Evol. 2022;54:20. et al. Deblur rapidly resolves single-nucleotide community sequence 23. Vernocchi P, del Chierico F, Putignani L. Gut microbiota profiling: patterns. mSystems. 2022;2:e00191-16. metabolomics based approach to unravel compounds affecting 46. Edgar RC, Flyvbjerg H. Error filtering, pair assembly and error cor‑ human health. Front Microbiol. 2016;7:1144. rection for next-generation sequencing reads. Bioinformatics. 24. Reverter A, Ballester M, Alexandre PA, Mármol-Sánchez E, Dalmau 2015;31:3476–82. A, Quintanilla R, et al. A gene co-association network regulating 47. Drake LE, Cuff JP, Young RE, Marchbank A, Chadwick EA, Symondson gut microbial communities in a Duroc pig population. Microbiome. WOC. An assessment of minimum sequence copy thresholds for iden‑ 2021;9:52. tifying and reducing the prevalence of artefacts in dietary metabarcod‑ 25. Ramayo-Caldas Y, Prenafeta-Boldú F, Zingaretti LM, Gonzalez-Rodriguez ing data. Methods Ecol Evol. 2022;13:694–710. O, Dalmau A, Quintanilla R, et al. Gut eukaryotic communities in pigs: 48. DeSantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, et al. diversity, composition and host genetics contribution. Anim Microbi‑ Greengenes, a chimera-checked 16S rRNA gene database and work‑ ome. 2020;2:18. bench compatible with ARB. Appl Environ Microbiol. 2006;72:5069–72. 26. Wang J, Chen L, Zhao N, Xu X, Xu Y, Zhu B. Of genes and microbes: 49. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy solving the intricacies in host genomes. Protein Cell. 2018;9:446–61. and high throughput. Nucleic Acids Res. 2004;32:1792–7. 27. Khan AA, Yurkovetskiy L, O’Grady K, Pickard JM, de Pooter R, Antonop‑ 50. Michener CD, Sokal RR. A quantitative approach to a problem of clas‑ oulos DA, et al. Polymorphic immune mechanisms regulate commen‑ sification. Evolution (N Y). 1957;11:130–62. sal repertoire. Cell Rep. 2019;29:541-550.e4. 51. Tamura K, Stecher G, Kumar S. MEGA11: molecular evolutionary genet‑ 28. Benson AK, Kelly SA, Legge R, Ma F, Low SJ, Kim J, et al. Individuality in ics analysis version 11. Mol Biol Evol. 2021;38:3022–7. gut microbiota composition is a complex polygenic trait shaped by 52. R Core Team. R: a language and environment for statistical computing. multiple environmental and host genetic factors. Proc Natl Acad Sci Vienna, Austria; 2020. https://​www.R-​proje​ct.​org/ Accessed 25 May USA. 2010;107:18933–8. 2020. 29. Goodrich JK, Davenport ER, Clark AG, Ley RE. The relationship between 53. Ward JH. Hierarchical grouping to optimize an objective function. J Am the human genome and microbiome comes into view. Annu Rev Stat Assoc. 1963;58:236–44. Genet. 2017;51:413–33. 54. Habier D, Fernando RL, Kizilkaya K, Garrick DJ. Extension of the Bayesian 30. Crespo-Piazuelo D, Migura-Garcia L, Estellé J, Criado-Mesas L, Revilla alphabet for genomic selection. BMC Bioinformatics. 2011;12:186. M, Castelló A, et al. Association between the pig genome and its gut 55. Pérez P, de los Campos G. Genome-wide regression and prediction microbiota composition. Sci Rep. 2019;9:8791. with the BGLR statistical package. Genetics. 2014;198:483–95. 31. Kundu P, Blacher E, Elinav E, Pettersson S. Our gut microbiome: the 56. Dietterich TG. Ensemble methods in machine learning. In: Multiple evolving inner self. Cell. 2017;171:1481–93. classifier systems. MCS 2000. Lecture notes in computer science, vol. 32. Gilbert JA, Quinn RA, Debelius J, Xu ZZ, Morton J, Garg N, et al. 18. Berlin: Springer; 2000. p. 1–15. Microbiome-wide association studies link dynamic microbial consortia 57. Azodi CB, Bolger E, McCarren A, Roantree M, de los Campos G, Shiu SH. to disease. Nature. 2016;535:94–103. Benchmarking parametric and machine learning models for genomic 33. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome prediction of complex traits. G3 (Bethesda). 2019;9:3691–702. datasets are compositional: and this is not optional. Front Microbiol. 58. Daetwyler HD, Pong-Wong R, Villanueva B, Woolliams JA. The impact of 2017;8:2224. genetic architecture on genome-wide evaluation methods. Genetics. 34. Khanal P, Maltecca C, Schwab C, Fix J, Bergamaschi M, Tiezzi F. Mod‑ 2010;185:1021–31. eling host–microbiome interactions for the prediction of meat quality 59. Lozupone C, Faust K, Raes J, Faith JJ, Frank DN, Zaneveld J, et al. Identi‑ and carcass composition traits in swine. Genet Sel Evol. 2020;52:41. fying genomic and metabolic features that can underlie early succes‑ 35. Buitenhuis B, Lassen J, Noel SJ, Plichta DR, Sørensen P, Difford GF, et al. sional and opportunistic lifestyles of human gut symbionts. Genome Impact of the rumen microbiome on milk fatty acid composition of Res. 2012;22:1974–84. Holstein cattle. Genet Sel Evol. 2019;51:23. 60. Faust K, Sathirapongsasuti JF, Izard J, Segata N, Gevers D, Raes J, et al. 36. Rothschild D, Weissbrod O, Barkan E, Kurilshikov A, Korem T, Zeevi D, Microbial co-occurrence relationships in the human microbiome. PLoS et al. Environment dominates over host genetics in shaping human Comput Biol. 2012;8: e1002606. gut microbiota. Nature. 2018;555:210–5. Calle‑García et al. Genetics Selection Evolution (2023) 55:29 Page 12 of 12 61. Bergamaschi M, Maltecca C, Schillebeeckx C, McNulty NP, Schwab C, Shull C, et al. Heritability and genome-wide association of swine gut microbiome features with growth and fatness parameters. Sci Rep. 2020;10:10134. 62. Xie H, Guo R, Zhong H, Feng Q, Lan Z, Qin B, et al. Shotgun metagen‑ omics of 250 adult twins reveals genetic and environmental impacts on the gut microbiome. Cell Syst. 2016;3:572-584.e3. 63. Saborío-Montero A, Gutiérrez-Rivas M, García-Rodríguez A, Atxaeran‑ dio R, Goiri I, López de Maturana E, et al. Structural equation models to disentangle the biological relationship between microbiota and complex traits: methane production in dairy cattle as a case of study. J Anim Breed Genet. 2020;137:36–48. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in pub‑ lished maps and institutional affiliations. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions