- Research Article
- Open Access

# Dominance and epistatic genetic variances for litter size in pigs using genomic models

- Zulma G. Vitezica
^{1}Email authorView ORCID ID profile, - Antonio Reverter
^{2}, - William Herring
^{3}and - Andres Legarra
^{4}

**Received:**19 April 2018**Accepted:**4 December 2018**Published:**22 December 2018

## Abstract

### Background

Epistatic genomic relationship matrices for interactions of any-order can be constructed using the Hadamard products of orthogonal additive and dominance genomic relationship matrices and standardization based on the trace of the resulting matrices. Variance components for litter size in pigs were estimated by Bayesian methods for five nested models with additive, dominance, and pairwise epistatic effects in a pig dataset, and including genomic inbreeding as a covariate.

### Results

Estimates of additive and non-additive (dominance and epistatic) variance components were obtained for litter size. The variance component estimates were empirically orthogonal, i.e. they did not change when fitting increasingly complex models. Most of the genetic variance was captured by non-epistatic effects, as expected. In the full model, estimates of dominance and total epistatic variances (additive-by-additive plus additive-by-dominance plus dominance-by-dominance), expressed as a proportion of the total phenotypic variance, were equal to 0.02 and 0.04, respectively. The estimate of broad-sense heritability for litter size (0.15) was almost twice that of the narrow-sense heritability (0.09). Ignoring inbreeding depression yielded upward biased estimates of dominance variance, while estimates of epistatic variances were only slightly affected.

### Conclusions

Epistatic variance components can be easily computed using genomic relationship matrices. Correct orthogonal definition of the relationship matrices resulted in orthogonal partition of genetic variance into additive, dominance, and epistatic components, but obtaining accurate variance component estimates remains an issue. Genomic models that include non-additive effects must also consider inbreeding depression in order to avoid upward bias of estimates of dominance variance. Inclusion of epistasis did not improve the accuracy of prediction of breeding values.

## Background

Genomics provides tools to understand the effects of genes and their interactions and new approaches for genetic improvement [1]. In quantitative genetics, partitioning genetic variance for a trait into statistical components due to additivity, dominance, and epistasis is useful for prediction and selection, even if it does not reflect the biological (or functional) effect of the underlying genes [2]. Additive, dominance and epistatic genetic variances of quantitative traits are required to estimate breeding values and for making optimal selection decisions.

Within-breed non-additive effects, and in particular epistasis, are often ignored in genetic improvement programs. However, the total genetic value of an animal is a function of both additive and non-additive effects and, taken together, these effects could result in better predictors of future phenotypes [3] and inform mate allocation. Indeed, assortative mating can improve the performance of livestock when dominance [4] and/or epistasis are/is present [5]. Evidence of non-additive variance in commercially important traits (e.g. for body depth in fish, [6]) has opened opportunities for specialized breeding schemes. To take the effects of dominance and/or epistasis into account, it is necessary to refine the estimation of non-additive variance components from phenotypic data.

The additive genetic effect is the part of an individual’s total genetic effect that is transmissible across generations from parents to offspring. In contrast, non-additive genetic effects (dominance and epistasis) can be seen as the residual genetic effect after fitting additive substitution effects. Although most of the genetic variance is additive [1, 7] and captures most of the functional epistatic action of genes, the epistatic variance should not be neglected. Knowing its magnitude in real data and exploring the predictive ability of a model that accounts for epistatic effects are of great relevance. It is also of interest to know how much gene-by-gene (G × G) interaction exists (e.g. to determine the effect of the same allele in different breeds) or to account for the contribution of epistasis to the creation of “new” additive variance over time [8].

Several authors have proposed the inclusion of dominance and epistatic effects in genetic evaluation using high-density marker genotypes (genomic evaluation) [5, 9–12]. Most epistatic models consider only additive-by-additive epistatic interactions (e.g. [9, 12]), although dominance-by-dominance and dominance-by-additive interactions may play a major role in heterosis and in inbreeding and outbreeding depression (see [13] p. 223). Recently, Vitezica et al. [14] proposed a flexible and general approach to construct “genomic” relationship matrices for populations that are or are not in Hardy–Weinberg equilibrium (HWE), e.g. F1 crosses. They proved that epistatic genomic relationship matrices for two or higher order interactions can be constructed using Hadamard products of additive and dominance genomic orthogonal relationships, regardless of the existence of HWE. They also pointed out that, nevertheless, standardization of genomic relationship matrices based on the trace of the relationship matrices is needed. However to date, models that use these relationship matrices have not been applied to real data.

Xiang et al. [15] proved analytically that, in the presence of directional dominance, inclusion of genomic inbreeding as a covariate is necessary to obtain correct estimates of dominance variance. Genomic inbreeding of an individual was shown to be correctly defined as the proportion of genotyped single nucleotide polymorphisms (SNPs) at which the individual is homozygous [16]. This was confirmed in real data by Xiang et al. [15] and Aliloo et al. [4]. The effect of fitting or not fitting genomic inbreeding on estimates of epistatic variance components is unknown.

In current pig production systems, the number of pigs weaned is a key factor to increase productivity, and litter size (e.g., total number of piglets born per litter) is one of the most important traits under selection, in maternal line breeding programs [17]. Although litter size in pigs has a low narrow-sense heritability [17–19], non-additive variance could still be abundant and should be ascertained. This is particularly important because the estimation of non-additive effects (dominance and epistasis) for litter size is of increasing interest since it can be used to define mate allocation strategies between selection candidates for developing new crossbreeding or even purebred breeding schemes.

The objectives of this work were to estimate additive and non-additive (dominance and epistasis) variance components for litter size for a real pig population, by considering or not inbreeding depression and to determine whether the accuracy of prediction of breeding values and total genetic values increases with the inclusion of dominance and epistatic effects.

## Methods

### Phenotypic and genomic data

Data for this study were provided by Genus plc (Hendersonville, TN, USA). Animal Care and Use Committee approval was not necessary for this study because the data were obtained from an existing database. Data on litter size (total number of piglets born per litter) were from a pig pure line. The average litter size was equal to 12.7 ± 3.1 and 13,369 records were available for 3619 sows. Genotypes for all sows were generated using the Illumina PorcineSNP60 BeadChip (Illumina, San Diego, CA). After quality control using default parameters by preGSf90 [20] (HWE, minor allele frequency, SNP call rate and animal call rate), 38,779 autosomal SNPs remained and were used to build genomic relationship matrices.

### Genomic evaluation models

In these analyses, we implicitly assume that allele substitution effects of quantitative trait loci (QTL) are distributed as \({\varvec{\upalpha}}{ \sim }N\left( {0,{\mathbf{I}}\sigma_{\alpha }^{2} } \right)\), dominance effects as \({\mathbf{d}}{ \sim }N\left( {0,{\mathbf{I}}\sigma_{d}^{2} } \right)\), additive-by-additive epistatic effects as \(\left( {{\varvec{\upalpha \upalpha }}} \right){ \sim }N\left( {0,{\mathbf{I}}\sigma_{{\left( {\alpha \alpha } \right)}}^{2} } \right)\), and similarly for additive-by-dominant \(\left( {{\mathbf{\alpha d}}} \right)\) and dominant-by-dominant \(\left( {{\mathbf{dd}}} \right)\) effects. Note that under the assumption of HWE, the GBLUP “breeding” model in Vitezica et al. [10] that accounts for additive and dominance effects is a particular case of the model to analyze epistasis in Vitezica et al. [14], which in turn is an extension of the NOIA (natural orthogonal interactions) QTL analysis model by Alvarez-Castro and Carlborg [23]. Note that a “statistical” or “breeding” model implies that the covariance between breeding values and dominance deviations is zero [24]. The NOIA model is orthogonal, which means that genetic effects are defined in such a way that the addition of other genetic effects (e.g. dominance) to a model does not change the definitions of genetic effects (e.g. additive) that were already included in the model [23]. For instance, the additive effect is always the regression of genotypic value on gene content; the dominance effect is a regression of the remainder on a measure of heterozygosity or identity at the genotype, and so on. The NOIA model of Alvarez-Castro and Carlborg [23] achieves orthogonality automatically due to its construction of the incidence matrices \({\mathbf{M}}\) and \({\mathbf{W}}\) and their Kronecker products; Vitezica et al. [14] proved that, for analyses at the individual level, these Kronecker products can be reformulated as Hadamard products.

Statistically, orthogonality means that inclusion of new terms in the model does not change estimates of already included effects in an infinitely large population. For instance, in practice, going from an additive to an additive plus dominant model should not change the estimates of additive variance components much. The advantage of using orthogonality in genetics and breeding is that it is the only way to carry out the estimation of breeding values (additive “statistical” effects) in an unambiguous manner i.e. such that they do not depend on other genetic terms that are fitted in the model.

Variance components were estimated for five nested models that added, in succession, additive effects (\(A\)), dominance effects (\(A + D\)), additive-by-additive genetic effects (\(A + D + AA\)), additive-by-dominance genetic effects (\(A + D + AA + AD\)), and dominance-by-dominance genetic effects (\(A + D + AA + AD + DD\)). Genetic variances (\(\sigma_{A}^{2}\), \(\sigma_{D}^{2}\), \(\sigma_{AA}^{2}\), \(\sigma_{AD}^{2}\) and \(\sigma_{DD}^{2}\)) were estimated by Bayesian methods using Gibbs sampling in the software gibbs2f90 [25], available at http://nce.ads.uga.edu/wiki/. In total, 200,000 iterations were run, discarding the first 10,000 and keeping every 100th sample. Convergence was checked by visual inspection of trace plots for the chains. Post-Gibbs analysis included estimation of the effective sample size and the deviance information criteria (DIC) as a goodness-of-fit criterion [26].

### Predictive ability

For comparison purposes, we also investigated the predicted ability of phenotypes of “new” sows for the three models as the correlation \(cor\left( {y^{*} ,\hat{y}} \right)\) [30] where \(y^{*}\) is the corrected phenotypic observation obtained from the “whole” dataset \(\left( {y^{*} = y - X\hat{\beta } - f\hat{b}} \right)\), and \(\hat{y}\) is the predicted corrected observation from the “partial” dataset, which is equal to the sum of the estimated genetic values \(\left( {\hat{g}} \right)\) and the estimated permanent environmental effect \(\left( {\widehat{pe}} \right)\). The estimated genetic value included the estimated additive genetic effects (\(A\) model), the sum of estimated additive and dominant genetic effects (\(\hat{g} = \widehat{{g_{A} }} + \widehat{{g_{D} }})\) in the \(A + D\) model, and the sum of all genetic effects \((\hat{g} = \widehat{{g_{A} }} + \widehat{{g_{D} }} + \widehat{{g_{AA} }} + \widehat{{g_{AD} }} + \widehat{{g_{DD} }}\)), in the \(A + D + AA + AD + DD\) model. The regression coefficients of \(y^{*}\) on \(\hat{y}\) were also estimated.

## Results and discussion

Estimates of epistatic variance ranged from 0.11 ± 0.11 to 0.14 ± 0.12 for the additive-by-additive component, from 0.11 ± 0.09 to 0.12 ± 0.09 for the additive-by-dominant component, and were equal to 0.09 ± 0.09 for the dominant-by-dominant component. As explained in many previous papers [1, 7], the magnitude of the epistatic variances is trivial compared to that of the additive variance. However, non-orthogonal models can result in exaggerated estimates of epistatic variances (e.g. Su et al. [9] and Muñoz et al. [31]). Estimates of epistatic variances had a large standard error in all models (\(A + D + AA\), \(A + D + AA + AD\) and \(A + D + AA + AD + AA\)), which illustrates the difficulties in obtaining good estimates of epistatic variances also from genomic information, even when there are only two-way interactions.

Estimates (and posterior standard deviation) of narrow sense heritability and of variance components for models that included genomic inbreeding and successively added additive effects (\(A\)), dominance effects (\(A + D\)), additive-by-additive effects (\(A + D + AA\)), additive-by-dominance effects (\(A + D + AA + AD\)), and dominance-by-dominance effects (\(A + D + AA + AD + DD\))

Model | \(h^{2}\) | \(d^{2}\)* | \(i^{2}\)** | \(\sigma_{pe}^{2}\) | \(\sigma_{e}^{2}\) |
---|---|---|---|---|---|

\(A\) | 0.095 (0.013) | 0.931 (0.110) | 7.049 (0.116) | ||

\(A + D\) | 0.093 (0.013) | 0.022 (0.013) | 0.755 (0.143) | 7.053 (0.117) | |

\(A + D + AA\) | 0.093 (0.013) | 0.020 (0.010) | 0.016 (0.015) | 0.633 (0.178) | 7.052 (0.117) |

\(A + D + AA + AD\) | 0.092 (0.013) | 0.020 (0.011) | 0.024 (0.014) | 0.572 (0.179) | 7.051 (0.118) |

\(A + D + AA + AD + DD\) | 0.092 (0.013) | 0.019 (0.012) | 0.038 (0.017) | 0.450 (0.184) | 7.054 (0.117) |

Inclusion of non-additive (dominance and epistasis) effects in the model did not have a large effect on estimates of residual variance, but they reduced the permanent environmental variance (Table 1), which shows that non-additive genetic effects are captured by the permanent environmental effects if they are not included explicitly in the model. Similar results were observed by Aliloo et al. [4] when dominance was included in an additive model.

Only second-order epistatic effects (e.g. additive-by-additive) were included in this study. Although it is tempting to fit high-order epistasis terms given the relative ease of computing the Hadamar products of relationship matrices, caution is needed. First, the products \({\mathbf{G}} \odot {\mathbf{G}} \odot {\mathbf{G}} \ldots\) quickly tend to the identity matrix, in which case there is no hope of distinguishing genetic components from residual effects. Second, partitioning genetic variance into additive and non-additive statistical components does not indicate the importance of additive versus non-additive gene actions [2].

The global fit of the models was studied using DIC. Models with smaller DIC exhibit a better fit after accounting for model complexity. Differences in DIC between models of less than 7 units are considered as irrelevant [35]. The values of DIC were 68,365, 68,370, 68,367, 68,370 and 68,375 for the \(A\), \(A + D\), \(A + D + AA\), \(A + D + AA + AD\) and \(A + D + AA + AD + AA\) models, respectively, i.e. they were similar across models, and models that included dominance and epistasis do not appear to fit the data better than the simplest model.

The predictive ability of the models for selection candidates was assessed using two approaches: the three statistics based on method R and by conventional cross-validation. The method R statistics of bias \(\left( {b_{0} } \right)\), slope \(\left( {b_{1} } \right)\), and correlation \(\left( \rho \right)\) were based on comparison of EBV obtained from the whole data (up to 2014) and model \(A\) with EBV obtained from the partial data and models \(A\), \(A + D\) and \(A + D + AA + AD + AA\). Similar results (not shown here) were obtained when using the \(A + D\) and \(A + D + AA + AD + DD\) models in the whole dataset. For the correct models, \(b_{0} = 0\) and \(b_{1} = 1\) are expected. The estimate of \(b_{0}\) was equal to \(0.019\sigma_{A}\) across the models and the estimate of \(b_{1}\) was equal to 1.01, 1.04 and 1.09 in models \(A\), \(A + D\) and \(A + D + AA + AD + DD\), respectively. These estimates near zero for \(b_{0}\) and near to 1 for \(b_{1}\), suggest that the model is empirically unbiased. The statistic \(\rho\) (accuracy), measured as the correlation between the EBV of selection candidates based on “whole” and “partial” data, was around 0.73 for all models.

Predictive ability of phenotype for different models

Model | G | GD | GDI |
---|---|---|---|

\(cor\left( {y^{*} ,\hat{y}} \right)\)
| 0.135 | 0.136 | 0.136 |

Regression coefficient | 0.896 | 0.922 | 0.978 |

## Conclusions

Using orthogonal relationship matrices, empirically orthogonal estimates of additive, dominance and epistatic variances were obtained for litter size in a pig dataset. The broad-sense heritability for litter size was almost twice the narrow-sense heritability. Genomic models that include non-additive effects must consider simultaneously inbreeding depression based on inbreeding in order to obtain unbiased estimates of variance components. Inclusion of epistasis or dominance did not improve the accuracy of prediction of breeding or genotypic values.

## Declarations

### Authors’ contributions

ZV and AL proposed the models and methods of analysis. WH provided the pig data. ZV conducted the statistical analyses and wrote the first draft of the manuscript. AR and WH discussed the results, made suggestions and corrections. All authors read and approved the final manuscript.

### Acknowledgements

The authors thank INRA and CSIRO people for hosting their respective visits and for the friendly environment. They also thank Genus plc. (Hendersonville, TN, USA) for the data.

### Competing interests

The authors declare that they have no competing interests.

### Availability of data and materials

The dataset used and analyzed in this study is available from the corresponding author on reasonable request and with the approval of Genus plc.

### Consent for publication

Not applicable.

### Ethics approval and consent to participate

Animal Care and Use Committee approval was not necessary for this study because the data were obtained from an existing database.

### Funding

This work was supported by the INRA SELGEN metaprogram (projects EpiSel and OptiMaGics) and an INRA-CSIRO linkage action. The project was partly supported by the Toulouse Midi-Pyrénées bioinformatics platform.

### Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access**This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- Mäki-Tanila A, Hill WG. Influence of gene interaction on complex trait variation with multilocus models. Genetics. 2014;198:355–67.View ArticleGoogle Scholar
- Huang W, Mackay TFC. The genetic architecture of quantitative traits cannot be inferred from variance component analysis. PLoS Genet. 2016;12:e1006421.View ArticleGoogle Scholar
- Aliloo H, Pryce JE, González-Recio O, Cocks BG, Hayes BJ. Accounting for dominance to improve genomic evaluations of dairy cows for fertility and milk production traits. Genet Sel Evol. 2016;48:8.View ArticleGoogle Scholar
- Aliloo H, Pryce JE, González-Recio O, Cocks BG, Goddard ME, Hayes BJ. Including nonadditive genetic effects in mating programs to maximize dairy farm profitability. J Dairy Sci. 2017;100:1203–22.View ArticleGoogle Scholar
- Toro MA, Varona L. A note on mate allocation for dominance handling in genomic selection. Genet Sel Evol. 2010;42:33.View ArticleGoogle Scholar
- Joshi R, Woolliams JA, Meuwissen THE, Gjøen HM. Maternal, dominance and additive genetic effects in Nile tilapia; influence on growth, fillet yield and body size traits. Heredity (Edinb). 2018;120:452–62.View ArticleGoogle Scholar
- Hill WG, Goddard ME, Visscher PM. Data and theory point to mainly additive genetic variance for complex traits. PLoS Genet. 2008;4:e1000008.View ArticleGoogle Scholar
- Hill WG. “Conversion” of epistatic into additive genetic variance in finite populations and possible impact on long-term selection response. J Anim Breed Genet. 2017;134:196–201.View ArticleGoogle Scholar
- Su G, Christensen OF, Ostersen T, Henryon M, Lund MS. Estimating additive and non-additive genetic variances and predicting genetic merits using genome-wide dense single nucleotide polymorphism markers. PLoS One. 2012;7:e45293.View ArticleGoogle Scholar
- Vitezica ZG, Varona L, Legarra A. On the additive and dominant variance and covariance of individuals within the genomic selection scope. Genetics. 2013;195:1223–30.View ArticleGoogle Scholar
- Nishio M, Satoh M. Including dominance effects in the genomic BLUP method for genomic evaluation. PLoS One. 2014;9:e85792.View ArticleGoogle Scholar
- Jiang Y, Reif JC. Modeling epistasis in genomic selection. Genetics. 2015;201:759–68.View ArticleGoogle Scholar
- Lynch M, Walsh B. Genetics and analysis of quantitative traits. New York: Sinauer associates Inc.; 1998.Google Scholar
- Vitezica ZG, Legarra A, Toro MA, Varona L. Orthogonal estimates of variances for additive, dominance, and epistatic effects in populations. Genetics. 2017;206:1297–307.View ArticleGoogle Scholar
- Xiang T, Christensen OF, Vitezica ZG, Legarra A. Genomic evaluation by including dominance effects and inbreeding depression for purebred and crossbred performance with an application in pigs. Genet Sel Evol. 2016;48:92.View ArticleGoogle Scholar
- Silió L, Rodríguez MC, Fernández A, Barragán C, Benítez R, Óvilo C, et al. Measuring inbreeding and inbreeding depression on pig growth from pedigree or SNP-derived metrics. J Anim Breed Genet. 2013;130:349–60.PubMedGoogle Scholar
- Nielsen B, Su G, Lund MS, Madsen P. Selection for increased number of piglets at d 5 after farrowing has increased litter size and reduced piglet mortality. J Anim Sci. 2013;91:2575–82.View ArticleGoogle Scholar
- Guo X, Christensen OF, Ostersen T, Wang Y, Lund MS, Su G. Improving genetic evaluation of litter size and piglet mortality for both genotyped and nongenotyped individuals using a single-step method. J Anim Sci. 2015;93:503–12.View ArticleGoogle Scholar
- Bidanel JP. Biology and genetics of reproduction. In: Rothschild MF, Ruvinsky A, editors. The genetics of the pig. 2nd ed. London: AB International; 1998. p. 313–43.Google Scholar
- Aguilar I, Misztal I, Tsuruta S, Legarra A, Wang H. PREGSF90–POSTGSF90: computational tools for the implementation of single-step genomic selection and genome-wide association with ungenotyped individuals in BLUPF90 programs. In: Proceedings of the 10th world congress of genetics applied to livestock production: 17–22 Aug 2014, Vancouver.Google Scholar
- Hill WG, Mäki-Tanila A. Expected influence of linkage disequilibrium on genetic variance caused by dominance and epistasis on quantitative traits. J Anim Breed Genet. 2015;132:176–86.View ArticleGoogle Scholar
- VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91:4414–23.View ArticleGoogle Scholar
- Alvarez-Castro JM, Carlborg O. A unified model for functional and statistical epistasis and its application in quantitative trait loci analysis. Genetics. 2007;176:1151–67.View ArticleGoogle Scholar
- Falconer DS, Mackay TFC. Introduction to quantitative genetics. New York: Longman; 1996.Google Scholar
- Misztal I, Tsuruta S, Strabel T, Auvray B, Druet T, Lee DH. BLUPF90 and related programs (BGF90). In: Proceedings of the 7th world congress of genetics applied to livestock production, 19–23 Aug 2002, Montpellier. Communication No 28-07.Google Scholar
- Spiegelhalter DJ, Best NG, Carlin BP, Van Der Linde A. Bayesian measures of model complexity and fit. J R Stat Soc Ser B Stat Methodol. 2002;64:583–639.View ArticleGoogle Scholar
- Reverter A, Golden BL, Bourdon RM, Brinks JS. Method R variance components procedure: application on the simple breeding value model. J Anim Sci. 1994;72:2247–53.View ArticleGoogle Scholar
- Reverter A, Golden BL, Bourdon RM, Brinks JS. Technical note: detection of bias in genetic predictions. J Anim Sci. 1994;72:34–7.View ArticleGoogle Scholar
- Legarra A, Reverter A. Can we frame and understand cross-validation results in animal breeding? In: Proceedings of the 22nd conference of the association for the advancement of animal breeding and genetics, 2–5 July 2017, Townsville. p. 2273–80.Google Scholar
- Legarra A, Robert-Granié C, Manfredi E, Elsen JM. Performance of genomic selection in mice. Genetics. 2008;180:611–8.View ArticleGoogle Scholar
- Muñoz PR, Resende MFR, Gezan SA, Resende MDV, de los Campos G, Kirst M, et al. Unraveling additive from nonadditive effects using genomic relationship matrices. Genetics. 2014;198:1759–68.View ArticleGoogle Scholar
- Varona L, Sorensen D, Thompson R. Analysis of litter size and average litter weight in pigs using a recursive model. Genetics. 2007;177:1791–9.View ArticleGoogle Scholar
- Misztal I, Varona L, Culbertson M, Bertrand JK, Mabry J, Lawlor TJ, et al. Studies on the value of incorporating the effect of dominance in genetic evaluations of dairy cattle, beef cattle and swine. Biotechnol Agron Soc Environ. 1998;2:227–33.Google Scholar
- De Boer I, Hoeschele I. Genetic evaluation methods for populations with dominance and in breeding. Theor Appl Genet. 1993;86:245–58.View ArticleGoogle Scholar
- Plummer M, Best N, Cowles K, Vines K. CODA: convergence diagnosis and output analysis for MCMC. R News. 2006;6:7–11.Google Scholar