Varona et al. Genet Sel Evol (2019) 51:78 https://doi.org/10.1186/s12711-019-0521-3 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 A multivariate analysis with direct additive and inbreeding depression load effects Luis Varona1* , Juan Altarriba1, Carlos Moreno1, María Martínez‑Castillero2 and Joaquim Casellas3 Abstract Background: Inbreeding is caused by mating between related individuals and its most common consequence is inbreeding depression. Several studies have detected heterogeneity in inbreeding depression among founder indi‑ viduals, and recently a procedure for predicting hidden inbreeding depression loads associated with founders and the Mendelian sampling of non-founders has been developed. The objectives of our study were to expand this model to predict the inbreeding loads for all individuals in the pedigree and to estimate the covariance between the inbreed‑ ing loads and the additive genetic effects for the trait of interest. We tested the proposed approach with simulated data and with two datasets of records on weaning weight from the Spanish Pirenaica and Rubia Gallega beef cattle breeds. Results: The posterior estimates of the variance components with the simulated datasets did not differ significantly from the simulation parameters. In addition, the correlation between the predicted and simulated inbreeding loads were always positive and ranged from 0.27 to 0.82. The beef cattle datasets comprised 35,126 and 75,194 records on weights between 170 and 250 days of age, and pedigrees of 308,836 and 384,434 individual-sire-dam entries for the Pirenaica and Rubia Gallega breeds, respectively. The posterior mean estimates of the variance of inbreeding depres‑ sion loads were 29,967.8 and 28,222.4 for the Pirenaica and Rubia Gallega breeds, respectively. They were larger than those of the additive variance (695.0 and 439.8 for Pirenaica and Rubia Gallega, respectively), because they should be understood as the variance of the inbreeding depression achieved by a fully inbred (100%) descendant. Therefore, the inbreeding loads have to be rescaled for smaller inbreeding coefficients. In addition, a strong negative correlation (− 0.43 ± 0.10) between additive effects and inbreeding loads was detected in the Pirenaica, but not in the Rubia Gal‑ lega breed. Conclusions: The results of the simulation study confirmed the ability of the proposed procedure to predict inbreed‑ ing depression loads for all individuals in the populations. Furthermore, the results obtained from the two real data‑ sets confirmed the variability in the inbreeding depression loads in both breeds and suggested a negative correlation of the inbreeding loads with the additive genetic effects in the Pirenaica breed. Background which results in a reduction in the phenotypic yield of fit- Inbreeding is caused by mating between related indi- ness-related traits [2, 3], has been widely observed in ani- viduals and is associated with changes in the mean and mals, plants, and humans [2, 4] and is the most common variance of quantitative traits [1]. Inbreeding depression, consequence of inbreeding. The genetic basis of inbreed- ing depression stems from the high degree of homozygo- sity in inbred individuals, which reveals the presence of *Correspondence: lvarona@unizar.es recessive alleles or losses in the advantage of over-domi- 1 Departamento de Anatomía Embriología y Genética Animal, nance at heterozygous loci [2, 3]. Instituto Agroalimentario de Aragón (IA2), Universidad de Zaragoza, 50013 Saragossa, Spain Genetic variation within a population implies that Full list of author information is available at the end of the article inbreeding depression can vary depending on the © The Author(s) 2019. 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. Varona et al. Genet Sel Evol (2019) 51:78 Page 2 of 12 genotype of the individuals whose alleles produce iden- Therefore, it follows that the distribution of the inbreed- tity-by-descent (IBD) states in their progeny. Variability in ing loads is multivariate Gaussian as: inbreeding depression has been confirmed within ances-   tral lineages of Drosophila [5] and in sire families in dairy i ∼ N 0, Aσi2 , [6] and beef cattle [7]. Therefore, it is possible to define a specific hidden individual inbreeding depression load [8] where A is the numerator relationship matrix [13], and σi2 (hereafter referred to as inbreeding load), which can be is the additive genetic variance of the inbreeding loads. considered a hereditary trait [5, 9] with a phenotype that is only expressed when inbreeding occurs in its offspring. Inbreeding decomposition Lacy et al. [10] proposed a founder decomposition of Traditionally, inbreeding is split into several units that inbreeding, which assigns a different inbreeding depres- are attributed to founding individuals [10] by tracing sion to each founder of the population; yet, the com- back the origin of the alleles that might be IBD in inbred mon ancestors of inbred individuals are not restricted individuals. Alternatively, inbreeding can be decom- to the founding generation. Caballero and Toro [11] and posed into the sources of the co-ancestry between the García-Cortés et al. [12] proposed a Mendelian decom- parents of each individual [11, 12], which might include position of the inbreeding coefficient that assigns partial the founders of the population and the Mendelian sam- inbreeding coefficients to the founders and to the Men- pling of the non-founders. The difference between the delian sampling of the non-founders. This decomposi- two approaches is illustrated by a simple pedigree of five tion was the basis for the development of a mixed-model individuals (Fig. 1). Individual 5 has an inbreeding of approach that allows prediction of individual inbreeding 0.375. The founder decomposition splits it between two loads [8]. However, this model assumes that the founding founders, and attributes 0.25 and 0.125 to individuals 1 and non-founding effects on inbreeding depression are and 2, respectively. In contrast, the Mendelian decom- distributed independently, and cannot be used to predict position attributes 0.1875 and 0.0625 to founders 1 and the inbreeding load of individuals that do not have inbred 2, respectively, and 0.125 to the Mendelian sampling of progeny. In this study, we propose an alternative param- non-founder 3. Note that, under the founder decomposi- eterization that can predict the inbreeding load for each tion, the genealogical inbreeding is the sum of the partial individual in the pedigree, and can provide estimates coefficients from the founders, while under the Mende- of the covariance between the inbreeding loads and the lian decomposition it is the sum of the partial coefficients additive genetic effects for the trait of interest. We tested from the founders and from the Mendelian sampling of the procedure with a simulation study and two large the non-founders. datasets of records on weaning weight in the Pirenaica and Rubia Gallega beef cattle breeds. Model Based on the Mendelian decomposition of inbreeding, Methods Casellas [8] proposed a linear model that explains the Theory phenotypic data ( y ) with two vectors of random effects: Inbreeding load the additive genetic (a ), and one extra effect (εi ) that is The recessive alleles that appear in homozygosity caused attributed to the inbreeding depression generated by the by IBD from a common ancestor in the paternal and inbreeding loads of their ancestors (the founders and the maternal lineages can cause inbreeding depression. Each Mendelian sampling of the non-founders). Thus, ancestor might have different sets of recessive alleles; thus, the genotype of each individual determines an y = 1′ µ + Za + εi + e, (2) inbreeding load that can be interpreted as the effects on where µ is the general mean, e is the vector of the residu- the trait performance of its inbred descendants. Inbreed- als, 1 is a vector of ones with the same length of y, Z is the ing load can be considered a heritable trait [5, 9] and, incidence matrix and T is the matrix that contains the from the perspective of the ancestors, it acts additively partial inbreeding coefficients from the Mendelian because the alternative alleles never interact in their decomposition and connects the phenotypic data of inbred progeny. inbred individuals and the inbreeding loads of their com- Thus, the polygenic inbreeding load (ii ) for the i th indi- mon ancestors fromthe paternal and maternal vidual can be decomposed [8] as follows:  lineages. Furthermore, a ∼ N 0, Aσa2 and e ∼ N 0, Iσe2 , where A   ii = is + id + εi , (1) is the numerator relationship matrix and σa2 and σe2 are their associated variance components. Casellas [8] pro- where is and id are the inbreeding loads for its sire and posed that the distribution of εi is εi ∼ N 0, Iσi2 with σi2  dam, respectively, and εi is its Mendelian sampling. being the variance component of the inbreeding loads. Varona et al. Genet Sel Evol (2019) 51:78 Page 3 of 12 Fig. 1 Founder and Mendelian decomposition of the inbreeding. Partial inbreeding coefficients for the fifth individual with the Founder and Mendelian decompositions of inbreeding However, this assumption is incorrect because the vari- σa2 σai   where G = , and σai is the covariance between ance that is attributed to the inbreeding loads of the σai σi2 founders must be larger than the variance that is attrib- additive effects and inbreeding loads. uted to the Mendelian sampling of the non-founders. Thus, the appropriate distribution of εi is εi ∼ N 0, Qσi2 ,   Example where Q is a diagonal matrix with a value 1 for the ele- Given the pedigree of five individuals (Fig. 1), we assume ments corresponding  to founder individuals, and a value that individuals 3, 4, and 5 have phenotypes 113, 87, and of 12 1 − FS +F 2 D for the elements corresponding to the 96, respectively. Applying the method proposed by García- non-founder individuals [11]. FS and FD are the genealog- Cortés et al. [12] a partial inbreeding coefficient of 0.250 ical inbreeding coefficient of the sire and dam of the indi- is obtained for the animal 4 generated by individual 1, and vidual, respectively. three partial inbreeding coefficients (0.1875, 0.0625, and In this study, we reparametrize the model in terms of 0.125) for animal 5, generated by individuals 1, 2, and 3, individual inbreeding loads ( i ) that, given their additive respectively. Therefore, the mixed-model equations for the nature, can be represented as follows: implementation of the model defined by Eq. (3) requires the following vector and matrices: i = Pi + εi ,   113 where P is a matrix with a diagonal of 0 s and 0.5 in the y =  87 , elements that link an individual with its sire and dam, 96 such that it establishes a recurrent relationship between the inbreeding loads of each individual with their par-   ents, as in Eq. (1). As a result, εi = (I − P)i and Model (2) 0 0 1 0 0 can be rewritten as follows: Z = 0 0 0 1 0 , 0 0 0 0 1 y = 1′ µ + Za + Ki + e, (3)   where K = T(I − P) and i ∼ N 0, Aσi . As a result, a 2 0 0 0 0 0   multiple trait model that includes the genetic covariance K = T(I − P) =  0.25 0 0 0 0 , between additive effects and inbreeding loads can then 0.125 0 0.125 0 0 be defined as:       0 0 0 0 0 a 0 i ∼N 0 ,G ⊗ A , (4) where T =  0.25 0 0 0 0 , 0.1875 0.0625 0.125 0 0 Varona et al. Genet Sel Evol (2019) 51:78 Page 4 of 12 1 0 0 0 0 with solutions:    0 1 0 0 0   and I − P =  −0.5  −0.5 1 0 0 .  µ̂ 104.88    −0.5 0 −0.5 1 0  â1     −2.64  0 0 −0.5 −0.5 1  â2   4.37   â     3   1.96  In addition, the numerator relationship matrix ( A )  â     4   −9.81   â5   and its inverse ( A−1) are as follows:    =  −6.24 .  î     1   −4.50  1 0 0.500 0.750 0.625    î2   −13.49       0 1 0.500 0.250 0.375   î3   −9.37      A =  0.500 0.500 1 0.750 0.875 ,      −5.41   0.750 0.250 0.750 1.250 1   î4  −7.07 0.625 0.375 0.875 1 1.375 î5 and Simulation 2 0.5 −0.5 −1 0   A pedigree of 30,000 individuals arranged in six discrete  0.5 1.5 −1 0 0  −1 generations of 5000 individuals (2500 sires and 2500 A =  −0.5 −1 3.071 −0.429 −1.143 ,    −1 0 −0.429 2.571 −1.143  dams) was simulated. Each generation was obtained from 0 0 −1.143 −1.143 2.286 the random mating of 20 sires and 500 dams that were chosen randomly from the previous generation. Thus, nei- which leads to the following mixed-model equations: ther selection nor purging were simulated. Once the pedi-  ′ gree was obtained, the direct additive genetic effects (a) 1′ Z 1′ K   11 µ̂ and the inbreeding loads ( i ) were simulated based on the  Z′ 1 Z′ Z + A−1 g11 σe2 Z′ K + A−1 g12 σe2  â  multivariate distribution from Eq. (4). The vector of the K′ 1 K′ Z + A−1 g21 σe2 K′ K + A−1 g22 σe2 î phenotypic records ( y ) was generated by Eq. (3). The par- ′   1y tial inbreeding coefficients required for the T matrix were =  Z′ y , derived from the procedure of García-Cortés et al. [12]. K′ y Two cases (i) and (ii) were simulated based on the fol- lowing parameters: where g11, g12, g22 are the elements of the inverse (G−1) of the covariance matrix between the additive and the i. µ = 100, σa2 = 100, σi2 = 10,000, σai = − 500, r(a, inbreeding load effects. Assuming σe2 = 10, σa2 = 2, i) = − 0.5 and σe2 = 100, σi2 = 1, and σai = −0.25, then ii. µ = 100, σa2 = 100, σi2 = 1000, σai = 0, r(a, i) = 0.0 and  g11 g12   2 −0.25 −1  0.516 0.129  σe2 = 100, G−1 = = = , g21 g22 −0.25 1 0.129 1.032 σai where r(a, i) = . σa σi and the equations become the following:   µ̂     296  3 0 0 1 1 1 0.375 0 0.125 0 0  â1   0 1.032 0.258 −0.258 −0.516 0 0.258 0.0645 −0.0645 −0.129 0    0      â2      0 0.258 0.774 −0.516 0 0 0.0645 0.1935 −0.129 0 0 â3   0           113   1  −0.258 −0.516 2.5849 −0.2211 −0.5897 −0.0645 −0.129 0.3962 −0.0553 −0.1474    â4     −0.516 −0.2211 −0.5897 −0.0553 −0.1474   1    0 2.327 0.121 0 0.332    87  â5       =  96   1 0 0 −0.5897 −0.5897 2.1794 0.125 0 −0.022 −0.1474 0.2949  ,     î1    33.75   0.375  0.258 0.0645 −0.0645 0.121 0.125 2.142 0.516 −0.500 −1.032 0        0 0.0645 0.1935 −0.129 0 0 −0.500 1.548 −1.032 0 0      0    î2      0.125 −0.0645 −0.129 0.3962 −0.0552 −0.022 1.0465 −1.032 3.185 −0.4423 −1.1794    12   î3        0 −0.129 0 −0.0553 0.3317 −0.1474 −1.032 0 −0.4423 2.6537 −1.1794  0      0 0 0 −0.1474 −0.1474 0.2948 0 0 −1.1794 −1.1794 2.3589   î4   0 î5 Varona et al. Genet Sel Evol (2019) 51:78 Page 5 of 12 The simulated variances of inbreeding Table 1 Distribution of partial inbreeding coefficients  2  loads σi were  2 greater than the additive variances σa as they reflect (number and percentage) based on their magnitude the inbreeding depression achieved by a completely in the Pirenaica and Rubia Gallega breeds inbred (100%) descendent. The percentage of phenotypic Partial inbreeding Pirenaica Rubia Gallega variance explained by the inbreeding depression loads < 10-4 10,996,924 (68.30%) 1,564,740 (30.80%) depends on the magnitude of the partial inbreeding coef- 10-4-10-3 3,967,842 (24.65%) 2,559,733 (50.38%) ficients. Therefore, in the first simulation case (i), the 10-3-10-2 1,028,578 (6.39%) 823,410 (16.21%) variance of inbreeding load for a partial inbreeding coef- 10-1-10-2 126,238 (0.78%) 124,155 (2.44%) ficient of 0.10 will be equal to the additive variance (100), > 10-1 9792 (0.06%) 8.419 (0.17%) whereas in the second (ii), it will be a tenth part (10). The simulated datasets were analyzed under Model (3). We used a Bayesian approach which included uniform priors for the mean (µ) and for the variance components were very low (Table 1); 99.3 and 93.0%, and 97.4 and (G) and σe2, and a multivariate Gaussian distribution for 81.2% were lower than 0.01 and 0.001 in the Pirenaica a and i . The analysis was implemented through a Gibbs and Rubia Gallega breeds, respectively. Sampler [14] with two chains of 550,000 iterations, after Once the partial inbreeding coefficients from the Men- discarding the first 50,000 iterations. delian decomposition of inbreeding were available, the phenotypic data were analyzed with the following model: Beef cattle data y = tc + Xb + Wp + Za + Ki + e, The dataset included 35,126 (Pirenaica) and 75,194 (Rubia Gallega) records on live weights (one per individual) where t is the vector of age at recording, b is the vector of between 170 and 250 days of age, which were provided systematic effects, which included sex and age of the dam by the two breeding associations and had been recorded (15 levels), p is the vector of herd-year-season effects since 1989 and 1991, respectively, until 2018. The aver- (Pirenaica, 6503 levels and Rubia Gallega, 5251 levels), age (± standard deviation) weights in the Pirenaica and and a , i , and e are the vectors of direct genetic, inbreed- the Rubia Gallega breeds were 258.4 kg (± 58.0 kg) and ing load, and residual effects, respectively. Furthermore, c 281.4 kg (± 46.7 kg), respectively. The genealogical infor- is the covariate with age at recording, and X, W , Z, and K mation dated back to 1920 (Pirenaica) and 1955 (Rubia are the appropriate incidence matrices. Gallega), and comprised 308,836 and 384,434 individ- Under a hierarchical Bayesian analysis, it was assumed ual-sire-dam records, respectively. The average depth of that prior distributions for the permanent environmen- the pedigree of the individuals for which there was phe- tal, additive, and inbreeding load effects followed a multi- notypic information was 5.9 and 3.9 generations for the variate Gaussian distribution: Pirenaica and the Rubia Gallega breeds, respectively. In   the Pirenaica breed, 248,471 (80.5%) individuals were p ∼ N 0, Iσp2 , inbred and had an average inbreeding (± standard devia- tion) coefficient of 0.030 (± 0.056). In the Rubia Gallega where σp2 is the variance of the permanent environmental breed, inbreeding affected 176,048 (45.8%) individuals, effects and and the average inbreeding coefficient (± standard devia-     tion) was 0.020 (± 0.048). The evolutions of the average u ∼N 0 ,G ⊗ A , inbreeding and the percentage of inbred individuals from i 0 1975 to 2017 are shown in Additional file 1: Figures S1 and Additional file 2: Figure S2 for Pirenaica and Rubia as in Model (4). Gallega, respectively. The prior distribution for systematic effects were uni- The partial inbreeding coefficients from the Mende- form between − M and M, being M a very large value. lian decomposition of inbreeding were derived from the The prior  −2distributions for σp2 was an inverted Chi- square χ (s, v) distribution with parameters s = 0 and  procedure of García-Cortés et al. [12], which generated 16,099,374 and 5,080,457 coefficients from 8721 and v = −2, and the prior distribution for G were inverted 3601 ancestors in the Pirenaica and Rubia Gallega breeds, Chi-square and inverted Wishart (IW (S, v)) distribution respectively. The average partial inbreeding coefficients with: (± standard deviation) were 5.7× ­10−4 (± 4.2×10−3) in  0 0  the Pirenaica and 1.5×10−3 (± 6.9×10−3) in the Rubia S= , 0 0 Gallega breed. Most of the partial inbreeding coefficients Varona et al. Genet Sel Evol (2019) 51:78 Page 6 of 12 and v = −3, that reduces to a uniform distribution. The standard deviation. Figure 3 shows the plots of the Gibbs sampler [14] was implemented with two chains of simulated direct genetic and inbreeding load effects 550,000 iterations, after discarding the first 50,000 itera- and their predictions (posterior mean estimates) for tions. Convergence was confirmed by visual inspection sires and for individuals without progeny. In addition, of the chains (see Additional file 3: Figures S3 and Addi- the correlations between the simulated and predicted tional file 4: Figure S4) and by the use of the CODA pack- effects for sires and individuals without progeny are age in the R software [15]. shown. The correlations between the direct genetic effects and their predictions ranged from 0.76 (individ- Results uals without progeny in case (i)) to 0.98 (sires in case Simulation study (ii)), and were always higher than their counterparts The posterior distributions of σa2 , σi2 , r(a, i) and σe2 for the inbreeding loads. Nevertheless, the correlations obtained from the simulation study are shown in Fig. 2. between simulated and predicted inbreeding loads The variance components used in the simulation were were always positive, and ranged from 0.27 (individuals within the posterior mean ± 2 times the posterior without progeny in case (ii)) to 0.82 (sires in case (i)). Fig. 2 Posterior distributions of the variance components in the simulation study. Posterior distributions of the additive variance σa2, inbreeding load variance σi2, genetic correlation r(a, i) and residual variance σe2 from the two simulation cases with the values used in the simulation (vertical red line) Varona et al. Genet Sel Evol (2019) 51:78 Page 7 of 12 Fig. 3 Relationship between simulated and predicted additive and inbreeding load effects in the simulation study. Simulated additive and inbreeding load effects (x axis) and their predictions (y axis) and correlation coefficients between them (r) for sires (a, b, e and f) and individuals without progeny (c, d, g and h). a–d Correspond to the first case of simulation and e–h to the second case Beef cattle datasets correlation between additive genetic and inbreeding load The posterior mean estimates of the variance compo- effects were − 0.43 in Pirenaica and − 0.04 in Rubia Gal- nents in the Pirenaica and Rubia Gallega breeds are in lega. The posterior mean estimates of the permanent Table 2. The posterior mean estimates of the additive environmental and residual variances were 1035.2 kg2 genetic variance and inbreeding load variance compo- and 483.1 kg2 in Pirenaica, and 320.0 kg2 and 1018.2 kg2 nents were equal to 695.0 kg2 and 29,966.8 kg2 in Pire- in Rubia Gallega. naica, and 439.8 kg2 and 28,222.4 kg2 in Rubia Gallega, Figure 4 presents, for each breed, a bivariate plot of respectively. The posterior mean estimates of the genetic the predictions (posterior mean estimates) of the direct Varona et al. Genet Sel Evol (2019) 51:78 Page 8 of 12 Table 2 Posterior mean estimates (and posterior standard additive and inbreeding load effects and a histogram that deviation) of variance components in the Pirenaica represents the distribution of the predicted inbreeding and Rubia Gallega breeds loads. In the Pirenaica breed, the bivariate plot indicated Population a clear negative relationship between direct genetic and inbreeding load effects, which was not apparent in the Pirenaica Rubia Gallega Rubia Gallega breed. The average (and standard devia- σa2 695.016 (25.688) 439.803 (18.121) tion) of the predicted inbreeding loads were -69.3 (58.6) σi2 29,966.800 (5868.275) 28,222.360 (5454.273) kg and − 32.3 (41.7) kg in Pirenaica and Rubia Gallega r(a, i) − 0.429 (0.102) − 0.043 (0.087) breeds, respectively, and ranged from − 369.1 to 234.7 σp2 1035.209 (27.714) 320.023 (11.273) in Pirenaica and from − 374.1 to 266.1 in Rubia Gallega. σe2 483.060 (15.010) 1018.207 (12.612) The proportion of individuals in the Pirenaica and Rubia Gallega breeds that had a positive predicted inbreeding σa2 is the additive genetic variance, σi2 is the variance of individual inbreeding loads, r(a, i) is the correlation between additive genetic and individual load was equal to 10.5% and 22.5%, respectively. inbreeding loads, σp2 is the permanent environmental variance and σe2 is the residual variance Fig. 4 Relationship between predicted additive and inbreeding load effects and distribution of inbreeding load effects for weaning weight in Pirenaica and Rubia Gallega. Predictions of the additive (x axis) and inbreeding load (y axis) effects for weaning weight in Pirenaica (a) and Rubia Gallega (b) and histograms of predictions of inbreeding load effects for weaning weight in Pirenaica (c) and Rubia Gallega (d) Varona et al. Genet Sel Evol (2019) 51:78 Page 9 of 12 Discussion to inbreeding but have relatives that have done so. The Inbreeding load and inbreeding depression are heterogene- approach is illustrated with the results obtained from a ous among populations [5] or families [6, 7], as inbreeding simulation study, which used the same model of analysis depression is affected by the variants of recessive alleles in to generate the simulated datasets. In spite of the limited each population, or even in each ancestor that generates amount of data, the posterior distribution of the variance inbreeding. Variability in inbreeding depression has been components did not differ significantly from the simu- confirmed in several livestock populations [16–21] with the lated ones (Fig. 2), which suggests that there is enough approach proposed by Lacy et al. [10], which is consistent information to recover the variance of inbreeding loads with the results of previous studies [22–24] that revealed and their covariance with the direct additive effects. The heterogeneity in the contribution of genomic regions to most outstanding result of the simulation study is that it inbreeding depression in several traits and species. demonstrates the ability of prediction of the inbreeding Management of inbreeding to prevent the undesir- loads of individuals without progeny, since the correla- able consequences is an important aspect of conserva- tions between the predicted and simulated inbreeding tion genetics [25] and in livestock and plant breeding loads were 0.27 and 0.48 in the first (i) and second (ii) programs [26, 27]. Purging is a natural process of selec- case of simulation, respectively (Fig. 3). The higher cor- tion against recessive alleles that are expressed in the relation in the second case was caused by the larger σi2 homozygous state [28]. Purging during past generations used in the simulation. The correlations obtained from can be inferred in pedigreed populations [29–31], and the individuals without progeny were lower than those evidence of purging has been detected in natural [32–34] for the sires, which had many inbred progeny, and they and livestock [29, 35] populations. were also lower than the correlations between simulated Genetic variability of the inbreeding loads among indi- and predicted additive effects. Nonetheless, these results viduals within a population can be taken into account in might provide a basis for the development of “artificial” ‘artificial’ purging strategies, to avoid the use of those indi- purging strategies in order to avoid the undesirable con- viduals with poor inbreeding loads for reproductive pur- sequences of inbreeding depression. One approach might poses. However, it requires prediction of the inbreeding be to develop a selection index by weighting the breeding loads of young individuals without inbred progeny. The values and inbreeding loads in an appropriate way. approach proposed by Lacy et al. [10] is useful for describ- Implementation of our model involved the construction ing the heterogeneity of inbreeding depression among of the mixed-model equations, which followed the same founders, yet in most livestock populations, the founders rationale as that of the maternal animal model [36]. The live several generations before the current candidates for model attributes two genetic effects to each individual; selection and, therefore, the prediction of the inbreeding one, i.e. a that it is expressed directly in the individual phe- loads based on Lacy’s approach is not useful for practi- notype, and another, i.e. i that it is expressed only in the cal ‘artificial’ purging strategies. Casellas [8] proposed phenotype of its inbred descendants. In spite of the com- an approach that uses the Mendelian decomposition of plexity of the model, the results obtained from the simu- inbreeding and provides prediction of the inbreeding lation study (Fig. 2) demonstrated its ability to recover loads from both founders and the Mendelian sampling of the simulated parameters from the genealogical and phe- the non-founders. However, it can only provide prediction notypic information. Computationally, the additional of the inbreeding loads of the individuals that contribute complexities beyond the standard mixed-model are the to the inbreeding of individuals with known phenotypes. calculation of the partial inbreeding coefficients from the Moreover, its estimates of the variance of the inbreeding Mendelian decomposition of inbreeding and the calcula- loads are biased, since it ignores the reduction in variance tion and storage of the K′ K block. First, the tabular method due to the Mendelian segregation of the progeny. that was implemented for the inbreeding decomposition Although the approach proposed by Casellas [8] was described in [12] involves the construction of a square the basis of the approach proposed in our study, rather matrix with the dimension of the pedigree. Thus, for com- than using directly the partial inbreeding coefficients putational feasibility, the partial inbreeding coefficients derived from García-Cortes et al. [12], the coefficients of each individual were calculated sequentially from a were transformed by the I–P recursive transformation reduced pedigree that only included its ancestors. Second, in order to predict the total individual inbreeding load the K′ K matrix is more dense than is the Z′ Z block and it ( i ). The parametrization is crucial because it allows the requires more memory storage and computational time. assumption that individual inbreeding loads follow a Nevertheless, the contributions of the partial inbreeding multivariate Gaussian distribution and permits the use coefficients to the K′ K matrix were squared and, thus, the of the numerator relationship matrix ( A ) to predict the very small partial inbreeding coefficients generated almost inbreeding loads of individuals that have yet to contribute null values, which could be discarded. Varona et al. Genet Sel Evol (2019) 51:78 Page 10 of 12 The posterior distributions of the inbreeding load vari- to date there are no published studies that provide esti- ances for weights at 210 days in the Pirenaica and Rubia mates of the correlations between additive genetic and Gallega breeds indicated high variability among them, inbreeding load effects in livestock breeds, although Car- which is consistent with the variability of inbreeding olino and Gama [7] did not find a significant correlation depression among sire families reported by Carolino between the breeding values of the sires and the inbreed- and Gama [7]. The posterior mean estimates were much ing depression associated with them. The only study larger (29,966.8 and 28,222.4 in Pirenaica and Rubia Gal- available which provided an estimate of the genetic cor- lega breeds, respectively) than those of the other variance relation [40], although retracted, suggested a very small, components in the model, which was confirmed by the negative genetic correlation [41]. A possible explanation broad range of the prediction of the inbreeding loads. of this negative correlation is that older individuals may However, it should be noted that the model provides have a smaller additive genetic effect due to selection estimates of the inbreeding loads that must be under- and a higher inbreeding load as they are more exposed to stood as the inbreeding depression achieved by a fully purging than recent ancestors. This hypothesis may also inbred (100%) descendent. Logically, the predictions are explain the differences in the estimates of the genetic cor- not realistic and have to be rescaled. Thus, the expected relation between the Pirenaica and Rubia Gallega breeds, inbreeding depression from an ancestor that generates a since the depth of the pedigree in the Pirenaica was larger, partial inbreeding coefficient of 0.10 and has a predicted and, therefore, older individuals with more opportunities inbreeding load of − 100 kg is − 10 kg. of purging were included. Regardless of its cause, if the The analysis of the beef cattle datasets indicated negative correlation is confirmed, it will mean that indi- that, for most of the individuals, the prediction of their viduals with high breeding values tend to cause worse inbreeding load was negative, which is consistent with inbreeding depression if their descendants are inbred. several studies that have analyzed the inbreeding depres- Finally, it is worth noting that the assumed model uses sion in body weights of beef cattle [7, 37–39]. However, a a prior distribution of the inbreeding loads centered at significant proportion of the individuals (10.52% in Pire- zero. Further research is necessary to define alternative naica, and 22.54% in Rubia Gallega) had a positive predic- prior distributions that allow the mean of the inbreeding tion of their inbreeding loads. In spite of the prediction loads to differ from zero. error variance of the inbreeding loads, these figures may indicate that a relevant proportion of the individuals Conclusions might have had an incremental effect on the trait of inter- Our study proposes a mixed-model approach that est. Therefore, the potential inbreeding caused by mating includes individual inbreeding depression loads, direct between the descendants might not need to be avoided additive effects and the covariance between them. The or could even be favored, at least for the analyzed trait. It approach was applied to simulated data and two datasets should be noted that with a similar posterior mean esti- with records on weaning weight from beef cattle breeds mate of the variance component (29,966.8 and 28,222.4), (Pirenaica and Rubia Gallega) in Spain. The results the prediction of the inbreeding load effects was more demonstrated the ability of the model to recover the variable in the Pirenaica (SD: 58.6 kg) than in the Rubia simulated parameters and to provide prediction of the Gallega (SD: 41.7) breed. This was due to the different individual inbreeding loads of candidates for selection amounts of pedigree information available for the esti- without inbred progeny. mation of the inbreeding loads (5.9 vs. 3.9 generations in Pirenaica and Rubia Gallega breeds, respectively). A Supplementary information deeper pedigree generates more partial inbreeding coef- Supplementary information accompanies this paper at https​://doi. ficients (Table 1) and provides more reliable informa- org/10.1186/s1271​1-019-0521-3. tion for the prediction of the inbreeding loads, which is reflected in larger variances for the predictors. Additional file 1: Figure S1. Plot of the Gibbs Sampler chains for the inbreeding load variance in the Pirenaica Breed. The additive nature of the inbreeding loads implies that Additional file 2: Figure S2. Plot of the Gibbs Sampler chains for the it can be genetically correlated with other traits such as inbreeding load variance in the Rubia Gallega Breed. the direct additive effect on the analyzed traits. The mean Additional file 3: Figure S3. Evolution of the number of inbred and non- posterior estimate of the genetic correlation between the inbred individuals and the average percentage of inbreeding from 1975 to inbreeding loads and the direct additive effect in the Pire- 2017 in the Pirenaica Breed. naica breed was clearly negative (− 0.43) and the poste- Additional file 4: Figure S4. Evolution of the number of inbred and non- rior probability that it was negative was higher than 0.99. inbred individuals and the average percentage of inbreeding from 1975 to 2017 in the Rubia Gallega Breed. In contrast, the estimate of the genetic correlation in the Rubia Gallega breed was near zero. To our knowledge, Varona et al. Genet Sel Evol (2019) 51:78 Page 11 of 12 Acknowledgements 9. Moorad JA, Wade MJ. A genetic interpretation of the variation in inbreed‑ The authors thank the breeder associations of Pirenaica (Confederación ing depression. Genetics. 2005;170:1373–84. Nacional de Asociaciones de Ganado Pirenaico–CONASPI-) and Rubia Gallega 10. Lacy RC, Alaks G, Walsh A. Hierachical analysis of inbreeding depression in (Asociación Nacional de Criadores de Ganado Vacuno Selecto de Raza Rubia Peromyscus polionotus. Evolution. 1996;50:2187–200. Gallega—ACR​UGA​-). The authors also want to thank Bruce MacWhirter for his 11. Caballero A, Toro MA. Interrelations between effective population size help with language editing. and other pedigree tools for the management of conserved populations. Genet Res. 2000;75:331–43. Authors’ contributions 12. García-Cortés LA, Martínez-Ávila JC, Toro MA. Fine decomposition of the LV and JC proposed the model of analysis. LV and JC developed the FORTRAN90 inbreeding and the coancestry coefficients by using the tabular method. code. LV analyzed the data and wrote the first draft of the manuscript. JA and Conserv Genet. 2010;11:1945–52. CM provided the beef cattle data and created the databases. MMC and JC col‑ 13. Wright S. Coefficients of inbreeding and relationship. Am Nat. laborated in writing the final manuscript. All authors discussed the results, made 1922;56:330–8. suggestions and corrections. All authors read and approved the final manuscript. 14. Gelfand AE, Smith AFM. Sampling-based approaches to calculating marginal densities. J Am Stat Assoc. 1990;85:398–409. Funding 15. Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and This study was financed by two research contracts with the Confederación output analysis for MCMC. R News. 2006;6:7–11. Nacional de Asociaciones de Ganado Pirenaico–CONASPI- and the Asociación 16. Rodrigáñez J, Toro MA, Rodriguez MC, Silió L. Effect of founder allele Nacional de Criadores de Ganado Vacuno Selecto de Raza Rubia Gallega— survival and inbreeding depression on litter size in a closed line of Large ACR​UGA​-. White pigs. Anim Sci. 1998;67:573–82. 17. Casellas J, Piedrafita J, Caja G, Varona L. Analysis of founder-specific Availability of data and materials inbreeding depression on birth weight in Ripollesa lambs. J Anim Sci. The datasets and software used in the study are available from the corre‑ 2009;87:72–9. sponding author on reasonable request. 18. Casellas J, Vidal-Roqueta D, Flores E, Casellas-Vidal D, Llach-Vila M, Salgas-Fina R, et al. Epistasis for founder-specific inbreeding depression in Ethics approval and consent to participate rabbits. J Hered. 2011;102:157–64. Approval by an Animal Care and Use Committee was not required for this 19. Casellas J, Varona L, Ibáñez-Escriche N, Quintanilla R, Noguera JL. Skew study because the analyses were performed on pre-existing data that had distribution of founder-specific inbreeding depression effects on the been collected under standard farm management practices from commercial longevity of Landrace sows. Genet Res. 2008;90:499–508. breeders and administrated by the Pirenaica and Rubia Gallega breeding 20. Gulisija D, Gianola D, Weigel KA, Toro MA. Between-founder heterogene‑ associations within the context of their breeding schemes. ity in inbreeding depression for production in Jersey cows. Livest Sci. 2006;104:244–53. Consent for publication 21. Todd ET, Ho SYW, Thomson PC, Ang RA, Velie BD, Hamilton NA. Founder- Not applicable. specific inbreeding depression affects racing performance in Thorough‑ bred horses. Sci Rep. 2018;8:6167. Competing interests 22. Pryce JE, Haile-Mariam M, Goddard ME, Hayes BJ. Identification of The authors declare that they have no competing interests. genomic regions associated with inbreeding depression in Holstein and Jersey dairy cattle. Genet Sel Evol. 2014;46:71. Author details 23. Saura M, Fernández A, Varona L, Fernández AI, De Cara MÁR, Barragán C, 1 Departamento de Anatomía Embriología y Genética Animal, Instituto et al. Detecting inbreeding depression for reproductive traits in Iberian Agroalimentario de Aragón (IA2), Universidad de Zaragoza, 50013 Saragossa, pigs using genome-wide data. Genet Sel Evol. 2015;47:1. Spain. 2 Dipartimento di Agronomia Animali, Alimenti Risorce Naturali e 24. Ferenčaković M, Sölkner J, Kapš M, Curik I. Genome-wide mapping and Ambiente, Università degli Studi di Padova, 35122 Padua, Italy. 3 Departa‑ estimation of inbreeding depression of semen quality traits in a cattle ment de Ciència Animal i dels Aliments, Universitat Autònoma de Barcelona, population. J Dairy Sci. 2017;100:4721–30. 08193 Bellaterra, Spain. 25. Hedrick PW, Kalinowski ST. Inbreeding depression in conservation biol‑ ogy. Annu Rev Ecol Syst. 2000;31:139–62. Received: 8 March 2019 Accepted: 17 December 2019 26. Meuwissen TH. Maximizing the response of selection with a predefined rate of inbreeding. J Anim Sci. 1997;75:934–40. 27. Woolliams JA, Berg P, Dagnachew BS, Meuwissen THE. Genetic contribu‑ tions and their optimization. J Anim Breed Genet. 2015;132:89–99. 28. Crnokrak P, Barrett SCH. Perspective: purging the genetic load: a review of References the experimental evidence. Evolution. 2002;56:2347–58. 1. Falconer DS, Mackay TFC. Introduction to quantitative genetics. Harlow: 29. Gulisija D, Crow JF. Inferring purging from pedigree data. Evolution. Longman Group; 1996. 2007;61:1043–51. 2. Charlesworth D, Willis JH. The genetics of inbreeding depression. Nat Rev 30. García-Dorado A, Wang J, López-Cortegano E. Predictive model and Genet. 2009;10:783–96. software for inbreeding-purging analysis of pedigreed populations. G3 3. Hedrick PW, Garcia-Dorado A. Understanding inbreeding depression, (Bethesda). 2016;6:3593–601. purging, and genetic rescue. Trends Ecol Evol. 2016;31:940–52. 31. López-Cortegano E, Bersabé D, Wang J, García-Dorado A. Detection of 4. Leroy G. Inbreeding depression in livestock species: review and meta- genetic purging and predictive value of purging parameters estimated in analysis. Anim Genet. 2014;45:618–28. pedigreed populations. Heredity. 2018;121:38–51. 5. Fowler K, Whitlock MC. The variance in inbreeding depression and 32. Boakes EH, Wang J, Amos W. An investigation of inbreeding depression the recovery of fitness in bottlenecked populations. Proc Biol Sci. and purging in captive pedigreed populations. Heredity. 2007;98:172–82. 1999;266:2061–6. 33. Facon B, Hufbauer RA, Tayeh A, Loiseau A, Lombaert E, Vitalis R, et al. 6. Miglior F, Burnside EB, Hohenboken WD. Heterogeneity among families Inbreeding depression is purged in the invasive insect Harmonia axyridis. of Holstein cattle. In: Proceedings of the 5th world congress on genetics Curr Biol. 2011;21:424–7. applied to livestock production: 7–12 August 2014. Guelph; 1994. p. 34. Moreno E, Pérez-González J, Carranza J, Moya-Laraño J. Better fitness in 479–82. captive Cuvier’s Gazelle despite inbreeding increase: evidence of purg‑ 7. Carolino N, Gama LT. Inbreeding depression on beef cattle traits: esti‑ ing? PLoS One. 2015;10:e0145111. mates, linearity of effects and heterogeneity among sire-families. Genet 35. Mc Parland S, Kearney F, Berry DP. Purging of inbreeding depression Sel Evol. 2009;40:511–27. within the Irish Holstein-Friesian population. Genet Sel Evol. 2009;41:16. 8. Casellas J. On individual-specific prediction of hidden inbreeding depres‑ 36. Henderson CR. Applications of linear models in animal breeding. Guelph: sion load. J Anim Breed Genet. 2018;135:37–44. University of Guelph; 1984. Varona et al. Genet Sel Evol (2019) 51:78 Page 12 of 12 37. Burrow HM. The effects of inbreeding in beef cattle. Anim Breed Abstr. 41. Tomkins JL, Penrose MA, Greeff J, Lebas NR. Retraction. Science. 1993;61:737–50. 2011;333:1220. 38. Pereira RJ, Santana ML Jr, Ayres DR, Bignardi AB, Menezes GR, Silva LO, et al. Inbreeding depression in Zebu cattle traits. J Anim Breed Genet. 2016;133:523–33. Publisher’s Note 39. Reverter A, Porto-Neto LR, Fortes MRS, Kasarapu P, de Cara MAR, Burrow Springer Nature remains neutral with regard to jurisdictional claims in pub‑ HM, et al. Genomic inbreeding depression for climatic adaptation of lished maps and institutional affiliations. tropical beef cattle. J Anim Sci. 2017;95:3809–21. 40. Tomkins JL, Penrose MA, Greeff J, LeBas NR. Additive genetic breeding values correlate with the load of partially deleterious mutations. Science. 2010;328:892–4. 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