Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 https://doi.org/10.1186/s12711-021-00658-7 Ge n e t i c s Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access Opportunities and limits of combining microbiome and genome data for complex trait prediction Miguel Pérez‑Enciso1,2,4* , Laura M. Zingaretti2,4, Yuliaxis Ramayo‑Caldas3 and Gustavo de los Campos4 Abstract Background: Analysis and prediction of complex traits using microbiome data combined with host genomic information is a topic of utmost interest. However, numerous questions remain to be answered: how useful can the microbiome be for complex trait prediction? Are estimates of microbiability reliable? Can the underlying biological links between the host’s genome, microbiome, and phenome be recovered? Methods: Here, we address these issues by (i) developing a novel simulation strategy that uses real microbiome and genotype data as inputs, and (ii) using variance-component approaches (Bayesian Reproducing Kernel Hilbert Space (RKHS) and Bayesian variable selection methods (Bayes C)) to quantify the proportion of phenotypic variance explained by the genome and the microbiome. The proposed simulation approach can mimic genetic links between the microbiome and genotype data by a permutation procedure that retains the distributional properties of the data. Results: Using real genotype and rumen microbiota abundances from dairy cattle, simulation results suggest that microbiome data can significantly improve the accuracy of phenotype predictions, regardless of whether some microbiota abundances are under direct genetic control by the host or not. This improvement depends logically on the microbiome being stable over time. Overall, random-effects linear methods appear robust for variance compo‑ nents estimation, in spite of the typically highly leptokurtic distribution of microbiota abundances. The predictive performance of Bayes C was higher but more sensitive to the number of causative effects than RKHS. Accuracy with Bayes C depended, in part, on the number of microorganisms’ taxa that influence the phenotype. Conclusions: While we conclude that, overall, genome-microbiome-links can be characterized using variance component estimates, we are less optimistic about the possibility of identifying the causative host genetic effects that affect microbiota abundances, which would require much larger sample sizes than are typically available for genome- microbiome-phenome studies. The R code to replicate the analyses is in https://​github.​com/​migue​lpere​zenci​so/​ simub​iome. Background microbiome has been estimated to quantify its influence The relevance of microbial ecosystems associated with and has been named ‘microbiability’ (b2) [7], in symme- humans and animals for health and production is now try with the classical ‘heritability’ (h2) concept [8]. Previ- widely recognized, e.g., [1–6]. The fraction of pheno- ously, the term "hologenome" was coined to describe the typic variance for a given trait that is explained by the joint action of the genome and the microbiome on a phe- notype [9]. Numerous microbes are responsible for diseases, e.g. *Correspondence: miguel.perez@uab.es 1 ICREA, Passeig de Lluís Companys 23, 08010 Barcelona, Spain sepsis in humans, and they have been used for diag- Full list of author information is available at the end of the article noses for many years [10]. Yet, a consequence of the © The Author(s) 2021. 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. Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 2 of 20 microbiability being larger than zero is that the whole Large-scale studies in humans suggest a predominant microbiome can be used to predict complex pheno- role of the environment in shaping the gut microbiome types, regardless of whether it is a disease or a produc- [11]. However, regardless of the relative importance of tion trait. This is an important issue since the use of genetic and environmental factors in shaping the micro- microbiome data has the potential to alter how medi- biota, microbiome composition per se can have a predic- cal diagnoses in humans or management and breeding tive value. Yet, the use of microbiota for the prediction decisions in agricultural species are performed. of future phenotypes or disease outcomes requires some Several studies have demonstrated the potential value level of stability of the microbiome over time. In the case of microbiome data for the prediction of complex- of the gastrointestinal tract, microbiota colonization traits. For example, Rothschild et al. [11] showed that starts at birth, when vertical transmission occurs through microbiome data can be used to improve accuracy in the mother’s birth canal. Then, microbiota diversity and the prediction of obesity and many other phenotypes in richness tend to increase as the host ages, to stabilize at humans. Likewise, Lloyd-Price et al. [12] showed that adulthood [34, 35]. In ruminants, the microbial popula- microbiome data can predict future outbursts of bowel tions that inhabit the rumen appear progressively after disease in humans. Various studies have shown the birth and partially persist throughout life [36]. power of microbiome data to predict methane emission As noted, the genome-microbiome-phenome is a com- and feed efficiency in cattle [4, 13–15], feed efficiency plex system but quantifying the relationships between and carcass traits in pigs [16, 17] and in poultry [18]. host-genome, microbiota, and phenotypes is important In addition, microbiota data from the rhizosphere has for the effective use of microbiome data for prediction been used to predict various plant phenotypes (e.g., of complex traits. Overall, although there are many pub- crop yield and diseases) [19]. Simultaneously, since the lished reports, we still lack detailed guidelines on the groundbreaking study of Meuwissen et al. [20], predic- joint use of microbiome and genome information for tion of complex traits using genomic information has the prediction of complex traits and on the reliability of been embraced in both plant [21] and animal breed- parameter inferences. The number of genes that affect ing [22], as well as in human genetics [23]. Therefore, microorganism abundance and that can be confidently combining the host’s genome and microbiome informa- identified, and the number of microorganism taxa that tion is a natural next step to improve the prediction of can influence a given phenotype remain unknown. With complex traits, a topic that is currently receiving much this work, our aim was to contribute to this important attention [16, 24]. topic by focusing on three inter-related questions: It is also important to realize that the composition of the microbiome can be affected by the host’s genome. 1. How useful can the microbiome be for prediction of Wang et al. [25] argued that it is evolutionarily justi- complex traits? fied for the microbiome to be under partial host genetic 2. Are microbiability estimates reliable? control since a non-negligible fraction of the cells in 3. Can the underlying biological genome-microbi- an adult body is made up of microbes, especially in the ome-links be inferred at a system’s level? On a more gut. Beginning with the seminal work by Pomp’s team refined level, the question that we aimed to address [26], several studies have confirmed the relationship is whether microbiome groups (e.g., OTU or gen- between the host’s genotype and microbiome composi- era) with sizable causal effects on phenotypes can be tion, e.g., [25, 27, 28]. These microbiome genome-wide identified with the typical size of current microbiome association studies (MWAS) suggest that microbiome data sets? abundances can be treated as any other complex trait in humans or livestock [27]. For example, Crespo-Piazuelo We address these questions via a novel simulation et al. [29] and Ramayo-Caldas et al. [30, 31] identified strategy that uses real microbiome and genotype data as several quantitative trait loci (QTL) that modulate bacte- inputs and by proposing a variance-component approach rial and eukaryotic communities in the gut of pigs and in that, in the spirit of mediation analyses, quantifies the rumen. Although the ‘heritability’ of individual amplicon proportion of phenotypic variance explained by the sequence variants (ASV) or operational taxonomic units genome and the microbiome. Importantly, the approach (OTU) is typically low, considering the whole microbi- allows simulation of a partial control of the host’s genome ome simultaneously should increase power of MWAS on the microbiome. This is accomplished using a partial [32]. In addition, although microbiome heritabilities permutation approach that preserves the distribution of vary according to the taxa level considered, they usually the genome and the microbiome. For the analyses, we increase as we move up from quasi-species to genus or used Bayesian Reproducing Kernel Hilbert Space (RKHS family levels, e.g., [33] [37]) and Bayes C [20] approaches. RKHS is similar to Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 3 of 20 genomic best linear unbiased prediction (GBLUP) [38], A novel data‑driven strategy to simulate while Bayes C is a variable selection approach that can microbiome‑genome‑phenotype experiments account for the possibility that some or all the features Ample literature and software are available on the simu- available in the genome or the microbiome have no effect lation of ‘standard’ complex phenotypes, e.g., [41–44]. on the trait of interest. We investigated the three above- However, these algorithms are not suited for some of mentioned questions across diverse causal scenarios that the scenarios presented in Fig. 1. Two issues make the examined the links between host genomes and microbi- simulation of the scenarios shown in Fig. 1 challenging: omes, and their relations with phenotype for a complex (i) microbiome data follow zero-inflated highly leptokur- trait. tic multivariate distributions [45, 46] and it is not obvi- ous how to sample from these distributions conditional on genome data, as is required in the ‘Recursive’ and Methods ‘Indirect’ scenarios; and (ii) in the absence of large-scale Causative scenarios considered published—and public—datasets, it is difficult to obtain Because the exact nature of the links between the accurate estimates of key parameters, such as microbia- genome (G), microbiome (B), and phenotype (y) is largely bility, to use in the simulations. To circumvent, or at least unknown and will likely vary from case to case, we used to alleviate, these constraints, we used publicly available the six generic causal models (‘scenarios’) illustrated in real data [4, 13] for both G and B. Fig. 1 to shed light on the nature of the genome-micro- Simulation under the ‘Joint’ scenario is straightforward biome-phenome links. In the ‘Null’ scenario, there is no since it assumes that G and B act independently (see link between any of the data-layers; while this scenario is below). Simulation under the ‘Recursive’ and ‘Indirect’ unlikely, it serves as an ‘overall null hypothesis’ and it is scenarios is not that obvious because causative abun- useful to assess potential biases in parameter estimates. dances are under genetic control and a link must exist The ‘Genome’ scenario assumes that only the genome between G and B. We solved this by rearranging abun- affects the phenotype. In turn, only the microbiome has dances within individuals such that the desired correla- a direct effect on the phenotype in the ‘Microbiome’ and tion between abundance and individual’s genotypes was ‘Indirect’ scenarios. In contrast to the ‘Microbiome’ sce- attained (see the Algorithm in Box 1 and the R-code at nario, the ‘Indirect’ scenario allows for some of the causa- https:// ​ g ithub. ​ c om/ ​ m igue ​ l pere ​ z enci ​ s o/ ​ S imub ​ i ome/​ tive abundances to be controlled genetically, which is blob/​master/​sortC​or.R). This strategy has the impor- similar to a scenario in which a phenotype is directly con- tant advantage that the distribution of abundances is not trolled by gene expression levels, and where gene expres- changed compared to the observed one. Figure 2 reca- sion is in turn genetically controlled [39, 40]. The ‘Joint’ pitulates the simulation strategy. The R code to replicate scenario is the simplest configuration for a trait that is the analyses is available at https://​github.​com/​migue​lpere​ under the influence of both the genome and the micro- zenci​so/​simub​iome). biome. It assumes that the microbiome and the genome are independent and that their effects on the phenotype are also independent. This is the most widely assumed Simulation details scenario, implicitly, or explicitly, in the literature, e.g., [4, We simulated the joint influence of the genome and the 11, 16]. The ‘Recursive’ scenario is similar to the ‘Joint’ microbiome on a quantitative trait by adding their contri- scenario but it accounts for the possibility that some butions plus random noise: causative OTU may be under partial genetic control by NQTN N the host. Therefore, in this scenario, the genome has both  OTU direct and indirect (microbiome-mediated) effects on the yi = αj zij + ωk xik + εi , (1) j=1 k=1 phenotypes. It should be noted that the ‘Recursive’ model does not assume that the same SNPs have both direct and where yi is the i -th individual record, αj is the genetic indirect effects, or that all OTU abundances are under effect of j-th causal SNP (QTN), with j = 1 to NQTN , genetic control. which is the number of QTN; zij is the genotype of the The causal models depicted in Fig. 1 were used to simu- i -th individual for the j-th SNP coded as -1, 0 and 1 late genome-microbiome-phenotype data under differ- (strict additivity was assumed for all QTN); ωk is the lin- ent configurations regarding the number of causative loci ear effect of the log-transformed abundance of the k -th (QTN) and the number of OTU with effects on the phe- OTU ( xik ), with k = 1 to NOTU , which is the number of notypes, as well as the number of OTU that are affected abundances that influence the phenotype; and εi is a nor- by the host’s genome. Tables 1 and 2 summarize the sim- mally distributed residual. The effect of an OTU can be ulation models and parameter values chosen. interpreted as the expected change in phenotype per unit Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 4 of 20 a b c d Fig. 1 Representation of the scenarios evaluated. G: genome, typically comprises SNP genotype data; B: microbiome; y: phenotype of interest; arrows indicate causality. An arrow from G to y indicates that there is a subset of elements of G (causative SNPs) that influence y; an arrow from G to B indicates that there is a subset of G that influences a subset of abundances in B, which, in turn, may also influence y. An arrow departing from B indicates there is a subset of microbial abundances (the causative abundances) that influence y. The SNPs that affect B need not necessarily be the same as SNPs that affect y directly in the Recursive scenario. Note that B can contain one or more sets of abundances such as archaea and bacteria communities, or different time or site sampling points. Without loss in generality, we assume B is a single community Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 5 of 20 Table 1 Definition of the evaluated scenarios and of the chosen parameters Scenario Abbreviation N QTN N OTU N OTU(g) h2 b2 Null 0 – – – 0 0 Genome G 100 0 0 r2 0 G500 500 0 0 r2 0 Microbiome M 0 25 0 0 r2 Indirect I 0 25 25 0 r2 Joint J 100 25 0 r 2/2 r 2/2 J500 500 25 0 r 2/2 r 2/2 Recursive R 100 25 25 r 2/2 r 2/2 R500 500 25 25 r 2/2 r 2/2 NQTN : number of SNPs with a direct causal effect on the phenotype y; NOTU : number of OTU with a direct effect on y; NOTU(g): number of OTU with a direct effect on y that are genetically determined, i.e., they are a subset of NOTU ; h2 is heritability, b2 is microbiability, and r 2 = h2 + b2. For r 2, values of 0.50 and 0.25 were considered. Causative OTU and SNPs were randomly sampled Table 2 Scenarios used to evaluate sensitivity of predictive accuracy to the number of causative OTU Scenario Abbreviation N QTN N OTU N OTU(g) r2 h2 b2 Joint J10 100 10 0 0.50 0.25 0.25 J100 100 100 0 0.50 0.25 0.25 J250 100 250 0 0.50 0.25 0.25 Recursive R10 100 10 5 0.50 0.25 0.25 R100 100 100 50 0.50 0.25 0.25 R250 100 250 125 0.50 0.25 0.25 Symbols are the same as in Table 1 Causative OTU and SNPs were randomly sampled increase in the log-transformed abundance of the OTU. G B G B G B Since abundances are on the log scale, this is equivalent to multiplicative effects of abundance on phenotype. y Equation (1) is valid for all scenarios NQTN in Fig. 1, except that y y the term involving markers j=1 αj zij is removed in Null Genome Microbiome the ‘Microbiome’ NOTU and ‘Indirect’ scenarios, while the term k=1 ωk xik is removed in the ‘Genome’ scenario. G B G B G B For the ‘Indirect’ and ‘Recursive’ scenarios, variation in abundances ( x) that is explained by the genome must also y y y be modeled (Fig. 1). Again, we can resort to a linear model Indirect Joint Recursive in which the log-transformed abundance is treated as a Fig. 2 Simulation scheme for the Recursive scenario (Fig. 1). a Real standard complex phenotype: input data comprises p genotypes (G matrix) and k taxa abundances NQTN (k) (B matrix). SNPs in grey are neutral, those in red act directly on the  phenotype y, and those in yellow/orange influence some OTU xik = βjk zij + ǫi , (2) abundances (marked in magenta color in B matrix); abundances j=1 in blue are not genetically controlled. b Given simulated effects, a genotypic value for abundance is obtained via Eq. (2). To obtain where xik is the log-transformed abundance level of the the required heritability, abundances in magenta are reordered; to k -th OTU that is under partial genetic control for the i simplify visualization, high abundances (represented by a darker -th individual, βjk is the genetic effect of the j-th QTN color) are associated with genotype ‘1’. A single SNP is shown on k -th abundance, and zij is the genotype of the i -th as causative for both OTU but there is no limit in practice. c The phenotype is simulated by adding the genome and the microbiome individual for the j-th SNP. The sum is across the QTN contributions plus a residual. d The general causal diagram that affect abundance of the k -th OTU, j = 1 to NQTN (k). Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 6 of 20 Note that abundances xik in Eq. (2) are a subset of those genome and microbiome is generated in the ‘Recursive’ in Eq. (1). Other non-causative abundances may also be and ‘Indirect ‘scenarios, mediated by Eqs. (1) and (2). under genetic control but that is irrelevant for our pur- A drawback of this algorithm is that it locally breaks pose. Based on this model, phenotype under the ‘Recur- the covariance between abundances of different OTU. sive’ scenario was simulated via a two-step procedure by To alleviate this, we permuted all abundances that fell first simulating abundances ( x ) using Eq. (2), followed by within the same OTU cluster. We clustered abundances simulating phenotype using Eq. (1) given the abundances using the R function hclust(dist(.), method = "ward.D2") obtained. and cut the tree in K = 500 clusters. We chose K = 500 We used real genome and microbiome data as input because the first quartile of the intra-cluster average cor- for the simulation procedure. We downloaded the abun- relation was above the third quartile of the average cor- dance of 4018 OTU from dairy cattle rumen (N = 750 relation between random abundances, i.e., clusters were [4]). A pseudo-count equal to one was added to zero made up of highly correlated abundances compared to abundances and all abundances were then total-sum average. We also explored K = 200 but we found no dif- scaled and log-transformed, which resulted in much less ference in predictive accuracy. To verify that the shuf- leptokurtic and asymmetric distributions than original fling algorithm did not alter the structure of the data, we raw abundances. In Eqs. (1) and (2), xik represents log- show the results of the principal component analysis of transformed abundances. the original microbiome set and a few shuffled microbi- High-density array genotypes for 750 Holstein cows ome sets in Additional file 1: Figure S2. Causative OTU were downloaded from [13]. To prune SNPs and facilitate were sampled from different clusters. computation, 35% of all SNPs (i.e. 32,204) were retained based on a minimum allele frequency of 0.01 and a maxi- Box 1 Finding a permutation of vectors x and y such mum missing genotypes percentage of 1%. The few miss- that the correlation between permuted vectors ing values were simply imputed with the mean. is a predetermined value ρ. Causative OTU and SNPs were randomly sampled, irre- Take x , y, and ρ, where x and y are arbitrary uncorre- spective of their frequency. Thirty simulation replicates lated vectors in Rn and 0 ≤ ρ ≤ 1 is the desired corre- per scenario were simulated. Under the ‘Joint’ scenario, lation. The aim is to find a permutation of y such that which assumes independence between G and B, we can the correlation cor(x, ysort ) = ρ, approximately. The simply sample the list of causative SNPs and abundances, algorithm can be equally applied when x and/or y con- simulate their effects, and apply Eq. 1 to generate phe- sist of integer numbers and normality is not required notype values given the observed genotypes and abun- either. The performance of the algorithm improves as dances. In the case of ‘Recursive’ and ‘Indirect’ scenarios, n increases and when normality does hold. it is not so obvious because we need to sample abun- dances that are under genetic control and a link must 1.Sort the values of x and y in increasing or decreasing exist between G and B (Eq. 2). We solved this issue by order. The correlation cor xsort , ysort ∼ = 1. rearranging abundances of a given OTU between individ- 2.Generate a dummy variable z =  sort + 2ewhere e y uals such that the desired correlation between abundance values are sampled from e ∼ N 0, Sy2 1−ρ ρ2 , with Sy2 and individual’s genotypes is attained. This strategy has being the sample variance of y . The correlation the important advantage that the distribution of abun- NQTNO cor(xsort , z) ∼ ρ. dances is not changed. Suppose γik = j=1 βjk zij 3.Create an index vector iy , which indicates how ysort is the simulated genetic effect of the i -th individual for should be reordered according log-transformed abundance of the k -th OTU (Eq. (2))   to the order of z. This dummy index iy = order y [order(z)] contains the and that the desired heritability for that abundance is h2k . order of y when values are back-sorted according to The algorithm (see Box 1) is based on the simple observa- the order of z. tion that, given any two vectors x and y , the correlation 4.Reorder iy = iy [rank(x)] to match the index with is maximum (ρ ~ 1) when observations in both vectors positions ysort in the original vector x . This is needed are sorted and ρ is ~ zero when they are shuffled. There- since x remains unchanged fore, there must be some order ysort that fulfills, approxi-    and only y is permuted. 5.The correlation cor x, y iy ∼ = ρ. mately, the constraint cor(x, ysort ) = ρ . For our purpose, we need to rearrange the observed abundances xk such that the correlation between the rearranged xk and γk is Parameter fitting hk , i.e. the square root of heritability for abundance of the Little is known on the number of OTU that influence a k -th OTU. The detailed algorithm is provided in Box 1. given phenotype and on how many of those are partly As a result, with this algorithm a covariance between Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 7 of 20 inherited. For that reason, we chose some extreme but Data analysis ‘educated’ values for each of the five scenarios depicted in It was not evident which predictive algorithm would work Fig. 1. We considered r 2 = h2g + h2b, where h2g and h2b are best for the complex scenarios simulated here, although the heritability and microbiability, respectively; r 2 = 0.25 results from the literature show that no approach is opti- is grossly the value reported by Difford [4] with N = 750, mal for all cases. Here, we compared the Bayes C algorithm whereas values closer to r 2 = 0.50 were reported by Wal- [20] and Bayesian RKHS regression, which is equivalent lace et al. [13] for some farms. In general, increasing r 2 to GBLUP [38], to assess prediction performance and reli- attempts to mimic the effect of increasing sample size. ability of parameter estimates. Both approaches represent We assumed h2g = h2b for the ‘Joint’ and ‘Recursive’ sce- extreme parameterizations in terms of priors applied and narios, as approximately reported by Difford et al. [4] were implemented in the BGLR R package [52]. To assess and Camarinha-Silva et al. [16]. The number of QTN predictive accuracy, 75 (10%) phenotypes were randomly was somewhat arbitrary and set to either 100 or 500, but removed and predicted with the fitted model using the the specific number of loci should not affect the results remaining data. The correlation between observed and much. Barton et al. [47] showed theoretically that most predicted phenotypes was used as a measure of predictive properties of the infinitesimal model hold as the number accuracy. of QTN increases even modestly (N > 20). The generic linear model used was: Numerous empirical and theoretical works have shown that genetic effects of QTN on phenotype are not uni- y = Zg + Wb + i + e, (3) formly distributed and can be approximated by a gamma- where y is the vector of the simulated phenotypes, g is like distribution [48, 49]. Thus, here we sampled direct the vector of SNP effects, Z is a matrix of the observed genetic effects α ∼ Ŵ (shape = 0.2, scale = 5), as suggested genotypes for the 33 k SNPs, b is the vector of the effects by Caballero et al. [50] and also used by us previously of log-transformed abundance of OTU, W is a matrix [51]. with all nb = 4018 log-transformed abundances for the Much less is known on the number of causative OTU n = 750 individuals, i contains the interaction between g ( NOTU ), although we can presume that NOTU should be and b, and e is a vector of residuals. Prior to the analy- smaller than the number of QTN. For instance, Duvallet ses, phenotypes, abundances and genotypic values were et al. [36] found in a large meta-analysis that the human standardized to a mean of zero and a SD of 1. For Bayes- diseases studied were affected by, on average, 10 to 15 ian RKHS, variance–covariance structures were specified changes in abundances at the genus level. Here, we con- to be Var g = G, Var(b) = B, and Var(i) = G ◦ B, with  sidered NOTU = 25 (0.6% of all OTU), although we also G = ZZ′ /n, B = WW ′ /nb, and ◦ denotes the direct prod- evaluated NOTU = 10, 100, and 250 (Table 2). Similarly, uct between matrices. The term i is intended to capture for the ‘Recursive’ and ‘Indirect’ scenarios, we took the any variance due to the interaction between genome and extreme scenario where all causative OTU are geneti- microbiome, similar to the additive x additive epistatic cally determined, i.e., NOTU = NOTU (g ). The genetic variance–covariance structure usually being obtained effects β on abundances (Eq. (2)) were sampled from the from G ◦ G [53, 54]. For priors, we used the default values same distribution β ∼ Ŵ(shape = 0.2, scale = 5) as the in the BGLR software. direct genetic effects α. Much less is known about the In Bayes C, we did not include an interaction explic- distribution of the effects of abundances, ω , on pheno- itly but, instead, computed a covariance between esti- type (Eq. (1)). We took as proxy the distribution of esti- mates b and g, as detailed below in Eq. 3. As priors π for mates of regression coefficients of methane emission on the probability of SNPs or abundances to enter the  Bayes abundances reported by Difford et al. [4] in their sup- C model, we used π ∼ Beta p0 = 5, π0 = 0.001 , which  plementary information S4, which can be approximated has expectation π0 and variance π0 (1 − π0 )/p0 + 1). by ~ Ŵ(shape = 1.4, scale = 3.8). Additional file 1: Figure We also considered a much  more liberal flat prior for S3 compares the distributions of QTN and OTU effects π ∼ Beta p0 = 2, π0 = 0.01 , but we did not observe and their fit to the empirical data of Difford et al. [4]. This strong differences (see Additional file 1: Figure S1). Unlike model predicts that the variance of the effects of OTU GBLUP, ‘heritability’ is not explicitly defined in a Bayes C on phenotype is on average wider and larger than that of framework but, here, we used the proposal by [52] (https://​ QTN. Although, at this point, this is speculative, it is sen- github.​com/​gdlc/​BGLR-R/​blob/​master/​inst/​md/​herit​abili​ sible to assume that only a few taxa have a sizeable influ- ty.​md) to estimate heritability and microbiability. In short, ence on a phenotype such as methane emission. Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 8 of 20 at each iteration i of the MCMC, the algorithm samples the (‘Indirect’ and ‘Recursive’ scenarios). In total, 50k itera- effects of the SNPs and OTU: tions, including 500 burn-in iterations and thinning every 5, were run for both the RKHS or Bayes C chains; a plot u(i) = Zĝ (i) , of the variances against iteration number indicated that convergence was attained with this number of iterations v (i) = W b̂(i) , (see Additional file 1: Figure S2). where u(i) and v (i) are genome and microbiome effects Results at the i -th iteration for the set of individuals, respec- How useful can the microbiome be for prediction tively, and ĝ (i) and b̂(i) are the sampled effects of the  SNPs of complex traits? and OTU abundances; therefore, Var(u (i) )/Var y and This logically depends on how much phenotypic variance Var(v (i) )/Var y are the sampled heritability and micro-   is jointly explained by the genome (h2) and the micro- biability in the i -th iterate, from which posterior means biome (b2), but also on how efficiently methods capture were estimated by averaging over iterations. For Bayes the relationship between the microbiome and the phe- Cgb, we also computed the sampled absolute covariance notype, and on how stable the microbiome is. It should between u and v for each iteration i , i.e.: be noted that prediction accuracy is conditionally inde- Niter   i i  pendent of heritability of the microbiome itself, i.e., given Cov u, y  = i=1 cov u , v the observed abundances B and observed genotypes G,  (4)    Var(y)Niter it does not matter whether the biological processes that generate B are affected by G. In other words, for a con- To assess how likely it is to identify causative OTU in stant r 2 = h2 + b2, prediction should not be affected by Bayes C, we computed the probability of a given OTU to whether the ‘Joint’ or ‘Recursive’ scenarios hold. Impli- enter the model over MCMC iterations. We ran a GWAS cations for genetic improvement, however, could be of abundances ( xk , k = 1, NOTU ) on SNP genotypes (zj , dramatically different. Breeding schemes that target the j = 1, NSNP ) using the R function lm(xk ∼ zj ) and com- microbiome could be designed provided the ‘Recur- puted the association P-value of both causative QTN, i.e., sive’ scenario holds but make no sense under the ‘Joint’ those that were simulated to affect abundances, and of scenario. neutral SNPs. This was done in the ‘Recursive’ scenario We compared the predictive performance of the Bayes- only, in which we also computed the heritabilities of all ian RKHS (GBLUP-like approach) and Bayes C [20] abundance levels using the RKHS model. Weakly inform- approaches when both genome and microbiome data ative priors for variances were used in this case to mimic were included in the model (Rgb and Cgb), including an a REML-like estimator. interaction term between g and b (Rgbx), or only genome We ran RKHS and Bayes C with complete models, i.e., data (Rg, Cg), or only microbiome data (Rb, Cb). For including genome, microbiome and their interaction, details of all models, see Table 3. First, we verified that and with partial models that considered only microbi- the null model, i.e., when phenotypes were permuted ome or genome information (Table 3). This was done relative to genotypes and abundances, did not result in to study confounding and to determine whether part of false predictive accuracies (see Additional file 1: Figure the variance in microbiome abundances was captured by S3A). Figures 3 and 4 show the predictive accuracies for the genome when the microbiome was partly heritable the two r 2 values considered, 0.50 and 0.25, respectively. In the case of r 2 = 0.50, we also explored the influence of Table 3 Statistical models used to analyze the data varying the number of causative OTU (Table 2). Overall, Bayes C showed better performance than Method Abbreviation Effects fitted RKHS, but it was more sensitive to an increase in the g b g×b number of QTN or of causative OTU (Fig. 3c). Impor- tantly, adding a term for the interaction between g and Bayesian RKHS Rgbx x x x (GBLUP) b did not improve the predictive performance of RKHS, Rgb x x - even in the ‘Recursive’ scenario (R). The full models Rgb Rg x - - and Cgb were much better than the partial models (Rg, Rb - x - Rb, Cg, Cb) when both h2 and b2 were larger than zero, Bayes C Cgb x x - as expected, i.e., for the ‘Joint’ and ‘Recursive’ scenarios Cg x - - (compare Fig. 3a vs. b and Fig. 4a vs. b). In these sce- Cb - x - narios, using both sources of variation improved pre- Effects included can be g (SNP genotypes), b (OTU abundances) and their diction compared to using only genome or microbiome interaction g × b Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 9 of 20 Fig. 3 Predictive accuracy with r2 = 0.50, computed as correlation between predicted and observed phenotypes across causal scenarios (Fig. 1), for each of the RKHS (Bayes C) models. Rgb (Cgb) considers microbiome and genome data; Rgbx includes genome and microbiome data and their interaction; Rg (Cg) includes genome data only; Rb (Cb) includes microbiome data only. a Complete models; b Partial models, c Effect of varying the number of causative OTU. Details of scenarios are in Tables 1 and 2: G: Genome; M: Microbiome; I: Indirect; J: Joint; R: Recursive. G500, J500 and R500 means 500 causative SNPs; J10 and R10, 10 causative OTU; J100 and R100, 100 causative OTU; J250 and R250, 250 causative OTU. Results are an average of 30 replicates per case data, especially when using Bayes C. Importantly, the also on how this variance is split between genome and predictive accuracy was slightly lower for the ‘Joint’ microbiome. Although this likely occurs because of the and ‘Recursive’ scenarios than for the ‘Microbiome’ larger noise in the ‘Recursive’ or ‘Joint’ scenarios than or ‘Genome’ scenarios. This indicates that the predic- in the ‘Microbiome’ or ‘Genome’ scenarios, it also sug- tive accuracy does not only depend on total r 2 , but gests that our strategy of analysis may not be optimal Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 10 of 20 Fig. 4 Predictive accuracy with r2 = 0.25, computed as correlation between predicted and observed phenotypes across causal scenarios (Fig. 1), for each of the RKHS (Bayes C) models. Rgb (Cgb) considers microbiome and genome data; Rgbx includes genome and microbiome data and their interaction; Rg (Cg) includes genome data only; Rb (Cb) includes microbiome data only. a Complete models; b Partial models. Details of scenarios are in Tables 1 and 2: G: Genome scenario; M: Microbiome; I: Indirect; J: Joint; R: Recursive. G500, J500 and R500 means 500 causative SNPs. Results are an average of 30 replicates per case and that there is room to develop more efficient tools, source of noise. This effect was less pronounced as r 2 especially for the ‘Recursive’ scenario. It should be increased. noted that the variance of prediction was larger for the We observed that the predictions were better when ‘Recursive’ than for the ‘Joint’ scenario for r 2 = 0.25, only the microbiome influenced the phenotype than i.e., heritability of abundances may be an additional when the genome was the only source of variation, a Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 11 of 20 phenomenon also observed with real data [13, 16, 24]. In Figure S3B). However, the behavior of estimates obtained this simulation, this likely occurred because the number with RKHS was different, as variance ratios are a priori of causative effects and of input variables (SNPs vs. OTU) bound between 0 and 1: average estimates of h2 and b2 was smaller for the ‘Microbiome’ or ‘Indirect’ scenarios were small yet non-zero when simulated values were than for the ‘Genome’ scenario. In fact, we observed a zero. Comparing Fig. 5a vs. Fig. 5b for estimates of h2 consistent negative correlation between the number of (b2) in the ‘Microbiome’ (‘Genome’) scenarios, estimates causative OTU and the predictive accuracy for both the are clearly zero with Bayes C, as expected, but not with ‘Joint’ and ‘Recursive’ scenarios (Fig. 3c). The number RKHS. The effect of prior information is much stronger of QTN also adversely affected prediction performance for low r 2, resulting in RKHS estimates that should be but mainly with Bayes C, whereas RKHS was not largely zero to be more biased (Fig. 6a). affected (see especially scenario G vs. G500 for h2 = 0.5 in It is interesting to observe that the heritability and Fig. 3a). microbiability estimates obtained with RKHS were less Taken together, our results suggest that predictive biased than those obtained with the Bayes C algorithm, accuracy could be increased by ~ 50% when consider- except when the true parameter is zero. This was more ing microbiome in addition to genome data, provided apparent for the ‘Joint’ and ‘Recursive’ scenarios and the microbiability is of the same order as the heritability r 2 = 0.5, as observed in the comparison between Fig. 5a (Fig. 3). This is probably an upper, optimistic limit, since and Fig. 5b. With Bayes C, an overestimation of b2 is it will be difficult to have microbiome data collected evident, regardless of the simulated value of r 2 . For the homogeneously over time and in different locations. ‘Joint’ and ‘Recursive’ scenarios, the upward bias in the While individuals can be genotyped at birth, the microbi- estimate of b2 was accompanied by an underestima- ome during early life stages may not be representative of tion of h2 , which indicates that variance estimates were that at adult or later stages. For instance, Maltecca et al. confounded when using the Bayes Cgb model (Fig. 5b). [55] showed that early life microbiota is not a good proxy However, this bias decreased when the number of for carcass composition in pigs, whereas later life micro- causative OTU increased. For instance, the bias in the biota is more strongly associated. b2 estimate was ~ 40% when NOTU = 10 but reduced We observed a roughly two-fold increase in predictive to ~ 10% with NOTU = 250 (Fig. 5c). In contrast, esti- accuracy when heritability was doubled for the ‘Genome’, mates obtained with RKHS were remarkably robust to ‘Joint’ and ‘Recursive’ scenarios, and a 50% increase for varying NOTU (Fig. 5d vs. e). Therefore, it is likely that the ‘Microbiome’ and ‘Indirect’ scenarios (Fig. 3 vs. the presence of a few causative OTU, but of large effect, Fig. 4). combined with the presence of highly leptokurtic abun- dance distributions, may result in biased parameter estimates when using Bayes C. It should be noted, in Are microbiability estimates reliable? turn, that estimates obtained with RKHS were inflated Reliable parameter estimates are needed to optimize the when they were actually zero, i.e., when the model design of breeding schemes, management practices or was overparameterized. This should be considered microbiome-wide association studies (MWAS [56]). They when interpreting microbiability estimates from real are also needed to understand the biology that underlies data. For instance, Difford et al. [4] report estimates of the interaction between microbiome and complex phe- h2 = 0.21 and b2 = 0.13 (N = 750) and found that G and notypes. To date, microbiability has usually been esti- B are independent. Assuming the number of causative mated using ‘standard’ linear methods, e.g., [4, 11, 32], OTU is small compared to the number of SNPs with much as we have done here. Thus, it is of interest to know an effect on abundances (QTN), our simulation results how accurate these estimates are. suggest that the estimate of b2 reported by Difford et al. Figures 5 (r 2 = 0.50) and 6 (r 2 = 0.25) show estimates of [2] might be inflated. If this is true, the actual contribu- the variance components for each of the scenarios and tion by the microbiome might be too small to improve model analyses from Tables 1, 2, 3. Bayes Cgb allows us to prediction over that obtained from using only genome assess whether h2 and/or b2 differ from zero: the microbi- data. Although Difford et al. [4] focused on inference ability estimate was on average near zero when the data rather than on prediction, they reported that no bac- were simulated according to the ‘Genome’ scenario and teria genera were significantly associated with meth- the heritability estimate was zero when the ‘Indirect’ or ane emissions. Other authors have reported multiple ‘Microbiome’ scenarios hold, as expected (Figs. 5b and microbial associations with methane emissions, includ- 6b). Similarly, estimates of both h2 and b2 were near ing members of bacterial, archaeal, fungal, and proto- zero when the null scenario held (see Additional file 1: zoan communities, e.g., [13, 30, 57–59]. Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 12 of 20 Fig. 5 Parameter estimates for r 2 = 0.50. Estimates of heritability (h2), microbiability (b2), and the correlation between the genome and microbiome ( gb) for each of the RKHS (Bayes C) analysis models: model Rgbx includes microbiome and genome data and their interaction; Cgb includes genome and microbiome data; Rg (Cg) includes genome data only, and Rb (Cb) includes microbiome data only. Details of simulation scenarios are in Tables 1 and 2: G, Genome scenario; M, Microbiome; I, Indirect; J, Joint; R, Recursive; G500, J500 and R500 mean 500 causative SNPs; J10 and R10, 10 causative OTU; J100 and R100, 100 causative OTU; J250 and R250, 250 causative OTU. The horizontal dashed lines indicate simulated h2 or b2 parameter values (0.25, 0.5 depending on the scenario). Results are an average of 30 replicates For comparison, panels c and d in Figs. 5 and 6 show again indicates that some confounding between b2 and estimates that were obtained when only genome or h2 occurred . In general, bias was lower when r 2 was microbiome information was used. The most noticeable greater but did not vanish. outcome is that bias in estimates of b2 was somewhat reduced relative to that found with Bayes Cgb, which Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 13 of 20 Fig. 6 Parameter estimates for r 2 = 0.25. Estimates of heritability (h2), microbiability (b2), and the correlation between the genome and microbiome ( gb) for each of the RKHS (Bayes C) analysis models: model Rgbx includes microbiome and genome data, and their interaction; Cgb includes genome and microbiome data; Rg (Cg) includes genome data only, and Rb (Cb) includes microbiome data only). Details of simulation scenarios are in Table 1: G: Genome scenario; M: Microbiome; I: Indirect; J: Joint; R: Recursive. The horizontal dashed lines indicate simulated h2 or b2 parameter values (0.125, 0.25 depending on the scenario). Results are an average of 30 replicates Can the underlying biological scenario and causative OTU some bias—the contribution of microbiability to the be recovered? phenotypic variance. It also appears possible to distin- An important goal of many experiments is to dissect the guish which of the ‘Microbiome’ or ‘Genome’ scenarios biological basis of the interactions between the microbi- fit a dataset best. Similarly, it appears possible to assess ome and the genome, even if this is not strictly needed whether both G and B contribute to the phenotypic vari- for prediction. So far, our simulations suggest that stand- ance, i.e., whether the ‘Recursive’ or ‘Joint’ scenarios are ard statistical methods can be used to quantify—with plausible. A question, however, is whether it is possible to Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 14 of 20 distinguish between the ‘Joint’ and ‘Recursive’ scenarios, to be performed, i.e., NOTU × NSNP . To illustrate the i.e., whether the data can indicate which of the ‘Indirect’ caveats of GWAS on abundances, Fig. 7a shows the dis- or ‘Microbiome’ scenarios is more plausible, if either. tribution of -log10 P-values of neutral SNPs vs. SNPs Furthermore, can causative OTU be identified? These with an effect on abundances. Taking the 5% empirical are far more difficult questions to answer than assess- threshold of the neutral P-value distribution to declare ing prediction performance or estimating microbiability. an association, simulations suggest that the P-values When the variance component estimates obtained under of only ~ 3% of the causative SNPs will be above that the ‘Joint’ and ‘Recursive’ scenarios are compared (Figs. 5 threshold, i.e., approximately what is expected by and 6), they appear to be nearly identical for the same r 2. chance. These P-values depend of course on the actual The two scenarios differ in that at least some causative number of causative SNPs and on abundance heritabili- OTU abundances are under partial genetic control in the ties, but most of the evidence to date points to a weak ‘Recursive’ scenario. Thus, the ‘Recursive’ scenario should relationship between the genome and the microbiome result in a covariance between G and B. For the RKHS [27]. It will be very difficult to identify causative SNPs modeling, we studied whether this covariance could for abundance using GWAS information alone [11, 25]. be partly captured by adding an interaction factor g ◦ b Another question of interest is what proportion of the (Eq. 3) much as imperfect disequilibrium can generate OTU that affect the phenotype we can expect to dis- ‘phantom’ epistasis [60]. However, the interaction esti- cover. One option is to count the frequency with which mates were not found to be greater in the ‘Recursive’ than a given OTU enters the Bayes C model during sampling. in the ‘Joint’ or ‘Microbiome’ scenarios when no interac- Figure 7b shows the probability of including a causa- tion was simulated (Figs. 5a and 6a). As for Bayes C, we tive OTU in the Bayes C sampling chain, which ranged investigated whether the two scenarios could be distin- from ~ 5% (b2 = 0.125) to ~ 20% (b2 = 0.25). On average, guished by analyzing the covariance cov(u(i) , v (i) )/var y about 50 (b2 = 0.25) and 30% (b2 = 0.125) of the causative   (see Methods). Again, these estimates were close to zero OTU were among the 5% most frequently included OTU regardless of the simulated scenario (Figs. 5b and 6b). in the Bayes C chain. Nevertheless, since the number An alternative approach to infer whether the ‘Recur- of causative OTU was 25, the rate of false positives was sive’ scenario holds or not is to run a genome-wide high. We can conjecture that only a few causative OTU association study (GWAS) for each OTU abundance. are likely to be identified in medium-sized experiments, If we identify significant SNPs for OTU that are likely such as simulated here. to influence the phenotype y, we could conclude that An alternative approach is a microbiome-wide asso- the ‘Recursive’ scenario is plausible. Unfortunately, this ciation study (MWAS), i.e., to perform a linear regres- analysis can be doomed by the large number of tests sion of the phenotype on each of the OTU abundances Fig. 7 Power in genome- (GWAS) and microbiome- (MWAS) wide association studies. a Distribution of -log10 P-values of a GWAS of abundances. b Probability of inclusion in the Bayes Cgb model of causative OTU for the two levels of microbiability considered. c Power of identifying a causative OTU computed as the probability of exceeding the 95% threshold of the empirical distribution of P-values in an MWAS for the Recursive scenario Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 15 of 20 and then select the significant OTU as potential causa- be subject to additive genetic variance in the experiment tive OTU [4]. Figure 7c shows the average power, defined reported by Difford et al. [4] (Fig. 9). as the percentage of true causative OTU among the 5% most significant OTU. In the ‘Recursive’ scenario, power Discussion was ~ 15 and ~ 20% for b2 = 0.125 and 0.25, respectively. Given the uncertainty with the true scenario, a flexible Again, this is not too satisfactory since we expect a high simulation approach is proposed fraction of false positives. In this scenario, it is perhaps Figure 1 represents highly simplified relationships more useful to consider probabilities of inclusion in the between the genome, microbiome, and phenotype. These Bayes C chain rather than P-values from linear regres- scenarios cover a wide range of possible interactions that sion since the former are the result of a joint analysis of are not exhaustive. For example, we did not consider the all OTU and can be used directly for prediction. case of a phenotype having a causative effect on abun- Finally, we investigated the pattern of abundance her- dances, i.e. the ‘reversed’ Microbiome scenario. The sce- itabilities. Figure 8a shows the simulated heritabilities narios chosen are nevertheless important to interpret for the causative inherited OTU, which approximately empirical data and can help to identify limiting factors followed a gamma distribution, as well as the estimated for prediction of complex traits. Furthermore, provided heritabilities for the causative OTU in the ‘Recursive’ sce- that a good fit is found, they will help the design of nario. The two distributions were rather similar, although experiments that combine microbiome and genetic data. the estimates were slightly shrunk towards zero, a con- We chose combinations of parameters that represent sequence of using a REML-like prior. Of course, a prob- extreme case scenarios and found that the results were, lem with real data is that we do not know which OTU are qualitatively, robust to the choice of parameters such inherited, and which are not, and the true distribution as r 2. However, a parameter that can be relevant is the of OTU heritability estimates will be a mixture due to number of causative microbiome taxa, i.e., those with an heritable and not heritable abundances. Figure 8b illus- effect on the phenotype, which appears to affect the bias trates the distributions of heritability estimates of neutral of microbiability estimates when Bayes C is used (Fig. 5). (non-inherited) and causative inherited OTU. In Fig. 8b, In this study, we proposed a new simulation procedure we mixed 1.7 neutral OTU per causative OTU, which that addresses some important challenges. First, the algo- is arbitrary since we do not know the actual number of rithm avoids the need for actual phenotype simulation OTU under genetic control, but the resulting mixture is by using real genotype and abundance data. Although similar to the distribution of heritabilities observed by we concede that this procedure may limit the general- Difford et al. [4] (Fig. 8c). If the distributions in Fig. 8b ity of the study, e.g., in terms of data size or abundance were representative of the true state of nature, this would patterns, we believe that the advantages of using real suggest that about 1/(1 + 1.7) ~ 40% of rumen OTU could data are numerous, since no simulation procedure can Fig. 8 Simulated, estimated and observed abundance heritabilities. a Distributions of simulated and GBLUP estimated heritabilities of abundance for causative OTU in the Recursive scenario. b Distribution of GBLUP estimated heritabilities of abundance for neutral and causative OTU in the Recursive scenario. c Distribution of estimates of heritabilities of OTU abundance reported by Difford et al. [4] Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 16 of 20 Fig. 9 Comparison of multiplicative (log) and threshold ‘Microbiome’ (M) and ‘Joint’ (J) scenarios (r 2 = 0.5) using Bayes Cgb. a Predictive accuracy, computed as the correlation between predicted and observed phenotypes; and b Estimates of heritability (h2) and microbiability (b2). Results are the average of 30 replicates. Scenarios M and J as specified in Table 1; the log transformation results are shown for completeness and are the same as in Figs. 3 and 5. Data are average of 30 replicates accommodate all known and unknown subtleties of the and statistical properties tend to depend primarily on the highly dimensional distributions at hand. Second, we actual read depth rather than on the specific ecosystem developed a permutation procedure (see Box 1) that [46]. allows linking previously uncorrelated data to fit a target Here, we have presumed that the effects on abun- correlation value. Also, by only permuting OTU within dances are additive on the log scale. Similar models are clusters, we minimized disruption of the overall covari- widely used in a diversity of scenarios. For example, ance structure (see Additional file 1: Figure S2). It should multiplicative models are used to accommodate fit- be noted that an OTU or SNP does not need to be causa- ness effects in evolutionary genetics [64] or to deal with tive to be useful for prediction and our simulation strat- highly leptokurtic distributions such as for abundances egy allows existing correlations between abundances or of microorganisms or of gene expression levels, which among SNPs (i.e. linkage disequilibrium) to be main- the log transformation addresses. In addition to the tained to mimic this fact. log-transformation, a widely popular choice in genet- For this work, we used real rumen microbiome and ics is the threshold model [8], which assumes the pres- genomic data from Holstein dairy cattle. Methane emis- ence of a continuous liability (here abundances) with an sions and the rumen ecosystem make a fundamental effect value of ‘0’ below a given threshold and ‘1’ other- area of microbiome research, but we believe our results wise. This model has the advantage of not depending on can be applied to other environments and species as whether abundances are log-transformed or not and is well, including humans, provided that similar generative also biologically sound since it is conceivable that a mini- parameters hold. We considered a range of values for mum microorganism abundance is required to trigger a heritability and microbiability, from 0 to 0.50, which cov- particular effect. To test the robustness of the log-trans- ers most of the estimates reported in the literature across formation, we simulated phenotypes such that 25% of the species, e.g., [3, 13, 30, 61]. Human studies tend to use a causative abundance observations were above the thresh- larger number of SNPs than the high-density cattle gen- old (i.e., abundances were binarized) and the analysis was otyping data used here. Yet, ample evidence shows that performed on the log transformed abundances, as before, the number of SNPs obeys a law of diminishing returns using Bayes C. As expected, using a ‘wrong’ model for the and a negligible difference in predictive performance has analyses was detrimental to prediction but not dramati- been observed between full sequence and high-density cally (Fig. 8a). Compared to the multiplicative model, genotyping data [51, 62, 63]. Distributions of OTU abun- parameter estimates were affected downwards (Fig. 8b). dances are always highly leptokurtic, as observed here, On this basis, we conjecture that the fundamental results obtained herein should hold even if the relationship Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 17 of 20 between OTU abundance and phenotype is not strictly There is growing evidence of interactions between the multiplicative. microbiome and the host genome [1, 25, 68] but it is not clear which approach is optimal to statistically model There is room for methodological developments this phenomenon. The interaction of the genome and the Numerous studies have reported estimates of microbi- microbiome poses challenging modeling and prediction ability for economically important traits, e.g., [4, 16, 30, problems. A main aim of this paper was to assess how 59], but the reliability of these estimates is not known. ‘standard’ linear approaches behave under these complex Estimates may be affected by the estimation proce- scenarios. Although the number of possible interactions dure and there are numerous alternatives to estimate to consider can be huge when the number of SNPs and b2, including Bayes C [20], GBLUP [38], and Bayesian the number of OTU is large, interactions between fea- RKHS regression using either Bray–Curtis dissimilarities tures in two high-dimensional sets can be modeled in as relationship matrix [30] or with the variance–covari- a Gaussian context using co-variance functions. These ance from the log-transformed OTU as kinship matrix functions are the Hadamard product of set-specific [14, 30]. Our results (Figs. 5 and 6) indicate that estimates similarity matrices such as the Hadamard product of a of microbiability obtained using BayesC may be biased SNP-derived and an OTU-derived ‘relationship’ matrix. upwards, especially when b2 is higher than 0.25 and the Such an approach has already been used to model, e.g. number of causative OTU is small. However, we found interactions between SNPs or between SNPs and envi- that estimates of b2 derived using Bayes C were very close ronmental covariates, e.g., [54], and here we employed to zero in the null scenario (see Additional file 1: Figure this in the RKHS framework through the G ◦ B covari- S1B). Thus, models using priors from the Spike-Slab fam- ance matrix. Unfortunately, no improvement in predic- ily, which consider a priori the possibility of null effects, tion was attained and the model was not able to capture can be used to test whether heritability or microbiabil- the covariance between genome and microbiome in the ity is substantial. In contrast, estimates obtained using recursive scenario. None of the methods evaluated here Bayesian RKHS were slightly less biased and less sensitive were optimal for both inference and prediction. Thus, to the number of causative OTU. However, estimates of there is room for improvement, but this requires proper variance fractions in the order ~ 0.10 were obtained even theoretical development that is beyond the scope of this if the true variance was zero. As a result, small effects work. Possibilities include extending recursive models as may not be distinguishable from null effects. Ramayo- done, e.g., by Saborio-Montero et al. [59], or improved Caldas [30] reported higher microbiability estimates selection indices, as in Weishaar [69]. using Bray–Curtis based kernels with Bayesian RKHS A wide array of penalized linear methods exists for than those using the log-transformed covariance matrix inference and prediction, which differ in their priors [70]. of abundances. Other methods have been proposed to We observed that Bayes C performed better than Bayes- select variables in a context of compositional data, e.g., ian RKHS in terms of prediction. We used a non-uniform [65, 66]. The behavior of estimation methods for microbi- distribution for simulating the effects of genes and OTU, ability certainly merits further research. specifically a gamma distribution, as is sustained by One conclusion from this work is that it will be difficult empirical and theoretical investigations, at least for gene to distinguish between some of the underlying scenarios effects [48, 50, 71]. Thus, it can be expected that the Bayes or to identify causative OTU and SNPs, at least by using C prior fitted the simulated data better than RKHS, for standard linear models, as was done here. The distinc- which the prior was flat across SNPs and OTU. The well- tion between ‘Joint’ and ‘Recursive’ scenarios is of special known ‘no-free-lunch theorem’ [72] in computer sciences relevance for breeding. The latter assumes partial genetic states that no method is superior in performance across control of some causative OTU. Yet, we found that both all scenarios or for all tasks. In agreement with this, we scenarios resulted in very similar patterns in terms of found that variance component estimates from RKHS predictive performance and parameter estimates (Figs. 3, were both slightly less biased and much less sensitive 4, 5, 6). Perhaps, a more powerful approach would be to the number of causative effects (OTU or SNPs) than to use structural equation models (SEM), which allow estimates from Bayes C. Previous work also showed that the inclusion of a variable both as an independent and Bayes C variance components are highly sensitive to the a dependent variable. Saborio-Montero et al. [67] com- genetic architecture of the trait [73]. pared a linear bivariate (one OTU and the phenotype) model with a SEM but found few differences between Final remarks: on using microbiome for prediction models. One limitation of their approach is that one SEM of complex traits was fitted for each abundance, instead of fitting several The utility of microbiome for the prediction of complex abundances simultaneously. traits, e.g., in prospective studies or in breeding, depends Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 18 of 20 crucially on its stability in time and space. Stability will Supplementary Information likely be lower for rare OTU than for core members [13], The online version contains supplementary material available at https://​doi.​ and so the abundance of causative OTU will be relevant org/​10.​1186/​s12711-​021-​00658-7. in this context. For instance, although measures of gas- trointestinal microbiome abundances are known to be Additional file 1: Figure S1. Comparison between flat and informative priors. Posterior distributions of heritability (black line) and microbiability repeatable, they cannot be expected to remain stable (blue dashed line) in a single replicate of the ‘Joint’ scenario, r 2 = 0.25. throughout an individual’s entire life span. After wean- Numbers in the panel titles are predictive accuracies for each prior and ing and under standard management conditions, e.g., method. Overall, mildly informative priors resulted in similar predictive accuracies as for flat priors and more reasonable posterior distributions. constant diet and absence of antibiotic treatment, the Figure S2. Gibbs sampling values of heritability (black line), microbiability diversity of monogastric gut microbiota increases with (red line), and correlation between genome and microbiome effects age of the host until its composition remains stable [34] (green line) in a single replicate of the ‘Recursive’ scenario (r 2 = 0.5). The comparison is between flat and informative priors. Figure S3. Results with Rumen microbial communities are highly resilient and the null model using the Bayes C model. This figure shows the results with host-specific [74, 75] but also change in early life, with the null model using the Bayes C model when samples were permuted the transition towards a more stable and adult-like rumi- relative to genotypes and abundances. (A) Prediction accuracy, computed as correlation between predicted and observed phenotypes, for each of nal ecosystem occurring between weaning and one year the Bayes C analyses: Cgb includes microbiome and genome; Cg includes of age [76]. Therefore, for prediction purposes, we rec- genome data only, and Cb includes microbiome data only. (B) Estimates ommend the inclusion of microbial data obtained at least of heritability (h2) and microbiability (b2) for each of the Bayes C analyses. Data are the average of 30 replicates. Figure S4. Plot of the principal after weaning, preferably at adulthood. Compared to component analysis of the original abundance data (top left) and three genomic data, this certainly limits the use of microbiota simulated datasets under the ‘Recursive’ model. Each dot corresponds to for prediction in breeding schemes. a single individual and data are log-transformed. The permutation has a negligible influence on the data structure. Figure S5. Distributions of effects. (A) Observed and fitted (red line) distribution of abundance linear Conclusions regression effects on methane emissions reported by Difford et al. [4]. (B) Comparison of gamma distributions used for sampling genetic (α) and To conclude, this study suggests that microbiome data OTUs’ (ω) effects: α ∼ Ŵ (k = 0.2, θ = 5) and ω ∼ Ŵ (k = 1.4, θ = 3.8), can significantly improve the prediction of complex phe- plotted in red and black lines, respectively (Eq. 1). notypes, regardless of whether some abundances are under direct genetic control or not. However, for this Acknowledgements strategy to be successful, medium- to large-sized experi- This work was developed during the stay of MPE and LMZ at the Michigan ments are required and the microbiome should be rela- State University kindly funded by GDLC. We appreciate insightful editing com‑ ments from editor J. Dekkers and the reviewers. tively stable and available prior to phenotype collection. This limits the use of the microbiome for prediction in Authors’ contributions breeding schemes as compared to genome data, which MPE, GDLC and LMZ conceived the research. MPE and LMZ performed the research. All authors discussed the research. MPE wrote the manuscript with can be collected at birth and remains unchanged. Never- the help from the other authors. All authors read and approved the final theless, important potential applications remain, such as manuscript. predicting methane emission in cattle, obesity and feed Funding efficiency, disease predisposition, or crop production LMZ is supported by a Ph.D. grant from the Ministry of Economy and Science using the soil metagenome. Overall, we show that stand- (MINECO, Spain), MPE is funded by Ministry of Science and Innovation-State ard linear methods can be used, in spite of the highly Research Agency (AEI) grant PID2019-108829RB-I00, and through the “Severo Ochoa Programme for Centres of Excellence in R&D” SEV‐2015‐0533 and leptokurtic distributions observed in OTU abundances. CEX2019-000902-S. YRC is recipient of a Ramon y Cajal fellowship (RYC2019- Given the specific advantages of each of the algorithms 027244-I) from the Spanish Ministry of Science and Innovation. evaluated, there is room for specific theoretical devel- Availability of data and materials opments that combine benefits from both. Nonetheless, GitHub site https://​github.​com/​migue​lpere​zenci​so/​simub​iomecontains the we argue that new models should be based on a better code, documentation and data links. understanding of the relation between the microbiome and the phenotype. It seems important to quantify, even Declarations approximately, the number of taxa that affect the pheno- Competing interests type and to characterize the distribution of their effects, The authors declare that they have no competing interests. as it may affect the reliability of parameter estimates (Fig. 5). However, we are far less optimistic (e.g., Fig. 7) Ethics approval and consent to participate Not applicable. with regards to the identification of causative OTU, and of the putative QTN that affect relative abundances. Author details 1 ICREA, Passeig de Lluís Companys 23, 08010 Barcelona, Spain. 2 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB, 08193 Bel‑ laterra, Barcelona, Spain. 3 Animal Breeding and Genetics Program, Institute Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 19 of 20 for Research and Technology in Food and Agriculture (IRTA), Torre Marimon, 20. Meuwissen THE, Hayes BJ, Goddard ME. Prediction of total genetic value 08140 Caldes de Montbui, Barcelona, Spain. 4 Dept. of Epidemiology & Bio‑ using genome-wide dense marker maps. Genetics. 2001;157:1819–29. statistics, and Dept. of Statistics & Probability, Michigan State University, East 21. Crossa J, Pérez-Rodríguez P, Cuevas J, Montesinos-López O, Jarquín D, de Lansing, MI 48824, USA. Los CG, et al. Genomic selection in plant breeding: Methods, models, and perspectives. Trends Plant Sci. 2017;22:961–75. Received: 25 January 2021 Accepted: 20 July 2021 22. Meuwissen T, Hayes B, Goddard M. Accelerating improvement of live‑ stock with genomic selection. Annu Rev Anim Biosci. 2013;1:221–37. 23. de Los CG, Gianola D, Allison DB. Predicting genetic predisposition in humans: the promise of whole-genome markers. Nat Rev Genet. 2010;11:880–6. References 24. Tierney BT, He Y, Church GM, Segal E, Kostic AD, Patel CJ. The predictive 1. Ruff WE, Greiling TM, Kriegel MA. Host–microbiota interactions in power of the microbiome exceeds that of genome-wide association immune-mediated diseases. Nat Rev Microbiol. 2020;18:521–38. studies in the discrimination of complex human disease. bioRxiv. 2020. 2. Zhang Q, Difford G, Sahana G, Løvendahl P, Lassen J, Lund MS, et al. https://​doi.​org/​10.​1101/​2019.​12.​31.​891978. Bayesian modelling reveals host genetics associated with rumen 25. Wang J, Chen L, Zhao N, Xu X, Xu Y, Zhu B. Of genes and microbes: solv‑ microbiota jointly influence methane emission in dairy cows. ISME J. ing the intricacies in host genomes. Protein Cell. 2018;9:446–61. 2020;14:2019–33. 26. Benson AK, Kelly SA, Legge R, Ma F, Low SJ, Kim J, et al. Individuality in 3. Maltecca C, Bergamaschi M, Tiezzi F. The interaction between microbi‑ gut microbiota composition is a complex polygenic trait shaped by ome and pig efficiency: a review. J Anim Breed Genet. 2020;137:4–13. multiple environmental and host genetic factors. Proc Natl Acad Sci USA. 4. Difford GF, Plichta DR, Løvendahl P, Lassen J, Noel SJ, Højberg O, et al. 2010;107:18933–8. Host genetics and the rumen microbiome jointly associate with methane 27. Goodrich JK, Davenport ER, Clark AG, Ley RE. The relationship between emissions in dairy cows. PLoS Genet. 2018;14:e1007580. the human genome and microbiome comes into view. Annu Rev Genet. 5. Kundu P, Blacher E, Elinav E, Pettersson S. Our gut microbiome: the evolv‑ 2017;51:413–33. ing Inner self. Cell. 2017;171:1481–93. 28. Blekhman R, Goodrich JK, Huang K, Sun Q, Bukowski R, Bell JT, et al. Host 6. Zhao H. The human microbiome and genetic disease: towards the genetic variation impacts microbiome composition across human body integration of metagenomic and multi-omics data. Hum Genet. sites. Genome Biol. 2015;16:191. 2021;140:701–2. 29. Crespo-Piazuelo D, Migura-Garcia L, Estellé J, Criado-Mesas L, Revilla M, 7. Difford GF, Lassen J, Løvendahl P. Genes and microbes, the next step in Castelló A, et al. Association between the pig genome and its gut micro‑ dairy cattle breeding. In Proceedings of the 6­ 7th Annual Meeting of the biota composition. Sci Rep. 2019;9:8791. European Federation of Animal Science: 29 August- 2 September 2016; 30. Ramayo-Caldas Y, Zingaretti L, Popova M, Estellé J, Bernard A, Pons N, Belfast; 2016. et al. Identification of rumen microbial biomarkers linked to methane 8. Falconer DS, Mackay TFC. Introduction to quantitative genetics. 4th ed. emission in Holstein dairy cows. J Anim Breed Genet. 2020;137:49–59. Harlow: Pearson Education Limited; 1996. 31. Ramayo-Caldas Y, Prenafeta-Boldú F, Zingaretti LM, Gonzalez-Rodriguez 9. Zilber-Rosenberg I, Rosenberg E. Role of microorganisms in the evolution O, Dalmau A, Quintanilla R, et al. Gut eukaryotic communities in pigs: of animals and plants: the hologenome theory of evolution. FEMS Micro‑ diversity, composition and host genetics contribution. Anim Microbiome. biol Rev. 2008;32:723–35. 2020;2:18. 10. Yanagihara K, Kitagawa Y, Tomonaga M, Tsukasaki K, Kohno S, Seki M, et al. 32. Weissbrod O, Rothschild D, Barkan E, Segal E. Host genetics and microbi‑ Evaluation of pathogen detection from clinical samples by real-time poly‑ ome associations through the lens of genome wide association studies. merase chain reaction using a sepsis pathogen DNA detection kit. Crit Curr Opin Microbiol. 2018;44:9–19. Care. 2010;14:R159. 33. Goorich JK, Davenport ER, Clark AG, Ley RE. The relationship between 11. Rothschild D, Weissbrod O, Barkan E, Kurilshikov A, Korem T, Zeevi D, the human genome and microbiome comes into view. Annu Rev Genet. et al. Environment dominates over host genetics in shaping human gut 2017;51:413. microbiota. Nature. 2018;555:210–5. 34. Bäckhed F, Roswall J, Peng Y, Feng Q, Jia H, Kovatcheva-Datchary P, et al. 12. Lloyd-Price J, Arze C, Ananthakrishnan AN, Schirmer M, Avila-Pacheco J, Dynamics and stabilization of the human gut microbiome during the first Poon TW, et al. Multi-omics of the gut microbial ecosystem in inflamma‑ year of life. Cell Host Microbe. 2015;17:690–703. tory bowel diseases. Nature. 2019;569:655–62. 35. Jakobsson HE, Abrahamsson TR, Jenmalm MC, Harris K, Quince C, 13. Wallace JR, Sasson G, Garnsworthy PC, Tapio I, Gregson E, Bani P, et al. Jernberg C, et al. Decreased gut microbiota diversity, delayed Bacteroi‑ A heritable subset of the core rumen microbiome dictates dairy cow detes colonisation and reduced Th1 responses in infants delivered by productivity and emissions. Sci Adv. 2019;5:eaav8391. caesarean section. Gut. 2014;63:559–66. 14. Ross EM, Moate PJ, Marett LC, Cocks BG, Hayes BJ. Metagenomic predic‑ 36. Furman O, Shenhav L, Sasson G, Kokou F, Honig H, Jacoby S, et al. tions: from microbiome to complex health and environmental pheno‑ Stochasticity constrained by deterministic effects of diet and age drive types in humans and cattle. PLoS One. 2013;8:e73056. rumen microbiome assembly dynamics. Nat Commun. 2020;11:1904. 15. Li F, Li C, Chen Y, Liu J, Zhang C, Irving B, et al. Host genetics influence the 37. de Los CG, Gianola D, Rosa GJM. Reproducing kernel Hilbert spaces rumen microbiota and heritable rumen microbial features associate with regression: a general framework for genetic evaluation. J Anim Sci. feed efficiency in cattle. Microbiome. 2019;7:92. 2009;87:1883–7. 16. Camarinha-Silva A, Maushammer M, Wellmann R, Vital M, Preuss S, Benne‑ 38. VanRaden PM. Efficient methods to compute genomic predictions. J witz J. Host genome influence on gut microbial composition and micro‑ Dairy Sci. 2008;91:4414–23. bial prediction of complex traits in pigs. Genetics. 2017;206:1637–44. 39. Schadt EE, Lamb J, Yang X, Zhu J, Edwards S, Thakurta DG, et al. An inte‑ 17. Khanal P, Maltecca C, Schwab C, Fix J, Tiezzi F. Microbiability of meat grative genomics approach to infer causal associations between gene quality and carcass composition traits in swine. J Anim Breed Genet. expression and disease. Nat Genet. 2005;37:710–7. 2020;138:223–36. 40. Hormozdiari F, van de Bunt M, Segrè AV, Li X, Joo JWJ, Bilow M, et al. 18. Vollmar S, Wellmann R, Borda-Molina D, Rodehutscord M, Camarinha- Colocalization of GWAS and eQTL signals detects target genes. Am J Hum Silva A, Bennewitz J. The gut microbial architecture of efficiency traits in Genet. 2016;99:1245–60. the domestic poultry model species Japanese quail (Coturnix japonica) 41. Peng B, Kimmel M. simuPOP: A forward-time population genetics simula‑ assessed by mixed linear models. G3 (Bethesda). 2020;10:2553–62. tion environment. Bioinformatics. 2005;21:3686–7. 19. Pereyra MA, Creus CM. Modifying the rhizosphere of agricultural crops to 42. Messer PW. SLiM: simulating evolution with selection and linkage. Genet‑ improve yield and sustainability: Azospirillum as a model rhizotroph. In: ics. 2013;194:1037–9. Mehnaz S, editor. Rhizotrophs: plant growth promotion to bioremedia‑ 43. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide tion. Microorganisms for sustainability, vol. 2. Singapore: Springer; 2017. p. complex trait analysis. Am J Hum Genet. 2011;88:76–82. 15–37. Pérez‑Enciso et al. Genet Sel Evol (2021) 53:65 Page 20 of 20 44. Pérez-Enciso M, Ramírez-Ayala LC, Zingaretti LM. SeqBreed: A python tool 62. Ober U, Ayroles JF, Stone EA, Richards S, Zhu D, Gibbs RA, et al. Using to evaluate genomic prediction in complex scenarios. Genet Sel Evol. whole-genome sequence data to predict quantitative trait phenotypes 2020;7:52. in Drosophila melanogaster. PLoS Genet. 2012;8:e1002685. 45. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome 63. van Binsbergen R, Calus MPL, Bink MCAM, van Eeuwijk FA, Schrooten datasets are compositional: and this is not optional. Front Microbiol. C, Veerkamp RF. Genomic prediction using imputed whole-genome 2017;8:2224. sequence data in Holstein Friesian cattle. Genet Sel Evol. 2015;47:71. 46. Weiss S, Xu ZZ, Peddada S, Amir A, Bittinger K, Gonzalez A, et al. Normali‑ 64. Crow J, Kimura M. An introduction to population genetics theory. New zation and microbial differential abundance strategies depend upon data York: Harper and Row; 1970. characteristics. Microbiome. 2017;5:27. 65. Susin A, Wang Y, Cao K-AL, Calle ML. Variable selection in microbiome 47. Barton NH, Etheridge AM, Véber A. The infinitesimal model: definition, compositional data analysis. NAR Genom Bioinform. 2020;2:Iqaa029. derivation, and implications. Theor Popul Biol. 2017;118:50–73. 66. Srinivasan A, Xue L, Zhan X. Compositional knockoff filter for high-dimen‑ 48. Eyre-Walker A, Keightley PD. The distribution of fitness effects of new sional regression analysis of microbiome data. Biometrics. 2020. https://​ mutations. Nat Rev Genet. 2007;8:610–8. doi.​org/​10.​1111/​biom.​13336 (in press). 49. Barton NH, Keightley PD. Understanding quantitative genetic variation. 67. Saborío-Montero A, Gutiérrez-Rivas M, García-Rodríguez A, Atxaerandio Nat Rev Genet. 2002;3:11–21. R, Goiri I, López de Maturana E, et al. Structural equation models to 50. Caballero A, Tenesa A, Keightley PD. The nature of genetic variation for disentangle the biological relationship between microbiota and complex complex traits revealed by GWAS and regional heritability mapping traits: methane production in dairy cattle as a case of study. J Anim Breed analyses. Genetics. 2015;201:1601–13. Genet. 2020;137:36–48. 51. Pérez-Enciso M, Forneris N, de Los CG, Legarra A. Evaluating sequence- 68. Kolde R, Franzosa EA, Rahnavard G, Hall AB, Vlamakis H, Stevens C, et al. based genomic prediction with an efficient new simulator. Genetics. Host genetic variation and its microbiome interactions within the Human 2017;205:939–53. Microbiome Project. Genome Med. 2018;10:6. 52. Pérez P, de Los CG. Genome-wide regression and prediction with the 69. Weishaar R, Wellmann R, Camarinha-Silva A, Rodehutscord M, Bennewitz BGLR statistical package. Genetics. 2014;198:483–95. J. Selecting the hologenome to breed for an improved feed efficiency in 53. Martini JWR, Toledo FH, Crossa J. On the approximation of interaction pigs—A novel selection index. J Anim Breed Genet. 2020;137:14–22. effect models by Hadamard powers of the additive genomic relationship. 70. Gianola D. Priors in whole-genome regression: The Bayesian alphabet Theor Popul Biol. 2020;132:16–23. returns. Genetics. 2013;194:573–96. 54. Muñoz PR, Resende MFR, Gezan SA, Resende MDV, de Los CG, Kirst M, 71. Avila V, Pérez-Figueroa A, Caballero A, Hill WG, García-Dorado A, López- et al. Unraveling additive from nonadditive effects using genomic rela‑ Fanjul C. The action of stabilizing selection, mutation, and drift on tionship matrices. Genetics. 2014;198:1759–68. epistatic quantitative traits. Evolution. 2014;68:1974–87. 55. Maltecca C, Lu D, Schillebeeckx C, McNulty NP, Schwab C, Shull C, et al. 72. Wolpert DH, Macready WG. No free lunch theorems for optimization. IEEE Predicting growth and carcass traits in swine using microbiome data and Trans Evol Comput. 1997;1:67–82. machine learning algorithms. Sci Rep. 2019;9:6574. 73. Berger S, Pérez-Rodríguez P, Veturi Y, Simianer H, de Los CG. Effective‑ 56. Power RA, Parkhill J, De Oliveira T. Microbial genome-wide association ness of shrinkage and variable selection methods for the prediction of studies: lessons from human GWAS. Nat Rev Genet. 2017;18:41–50. complex human traits using data from distantly related individuals. Ann 57. Roehe R, Dewhurst RJ, Duthie CA, Rooke JA, McKain N, Ross DW, et al. Hum Genet. 2015;79:122–35. Bovine host genetic variation influences rumen microbial methane 74. Cole NA. Effects of animal-to-animal exchange of ruminal contents on production with best selection criterion for low methane emitting and the feed intake and ruminal characteristics of fed and fasted lambs. J efficiently feed converting hosts based on metagenomic gene abun‑ Anim Sci. 1991;69:1795–803. dance. PLoS Genet. 2016;12:e1005846. 75. Weimer PJ. Redundancy, resilience, and host specificity of the ruminal 58. Huws SA, Creevey CJ, Oyama LB, Mizrahi I, Denman SE, Popova M, et al. microbiota: implications for engineering improved ruminal fermenta‑ Addressing global ruminant agricultural challenges through understand‑ tions. Front Microbiol. 2015;6:296. ing the rumen microbiome: past, present, and future. Front Microbiol. 76. Dill-Mcfarland KA, Breaker JD, Suen G. Microbial succession in the 2018;9:2161. gastrointestinal tract of dairy cows from 2 weeks to first lactation. Sci Rep. 59. Wallace RJ, Rooke JA, McKain N, Duthie CA, Hyslop JJ, Ross DW, et al. The 2017;7:40864. rumen microbial metagenome associated with high methane produc‑ tion in cattle. BMC Genomics. 2015;16:839. 60. de los Campos G, Sorensen DA, Toro MA. Imperfect linkage disequilib‑ Publisher’s Note rium generates phantom epistasis (& perils of big data). G3 (Bethesda). Springer Nature remains neutral with regard to jurisdictional claims in pub‑ 2019;9:1429–36. lished maps and institutional affiliations. 61. Goodrich JK, Waters JL, Poole AC, Sutter JL, Koren O, Blekhman R, et al. Human genetics shape the gut microbiome. Cell. 2014;159:789–99. 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