Parent-offspring genotyped trios unravelling genomic regions with gametic and genotypic epistatic transmission bias on the cattle genome

Several biological mechanisms affecting the sperm and ova fertility and viability at developmental stages of the reproductive cycle resulted in observable transmission ratio distortion (i.e., deviation from Mendelian expectations). Gene-by-gene interactions (or epistasis) could also potentially cause specific transmission ratio distortion patterns at different loci as unfavorable allelic combinations are under-represented, exhibiting deviation from Mendelian proportions. Here, we aimed to detect pairs of loci with epistatic transmission ratio distortion using 283,817 parent-offspring genotyped trios (sire-dam-offspring) of Holstein cattle. Allelic and genotypic parameterization for epistatic transmission ratio distortion were developed and implemented to scan the whole genome. Different epistatic transmission ratio distortion patterns were observed. Using genotypic models, 7, 19 and 6 pairs of genomic regions were found with decisive evidence with additive-by-additive, additive-by-dominance/dominance-by-additive and dominance-by-dominance effects, respectively. Using the allelic transmission ratio distortion model, more insight was gained in understanding the penetrance of single-locus distortions, revealing 17 pairs of SNPs. Scanning for the depletion of individuals carrying pairs of homozygous genotypes for unlinked loci, revealed 56 pairs of SNPs with recessive epistatic transmission ratio distortion patterns. The maximum number of expected homozygous offspring, with none of them observed, was 23. Finally, in this study, we identified candidate genomic regions harboring epistatic interactions with potential biological implications in economically important traits, such as reproduction.

1 Introduction

Transmission bias phenomenon (or transmission ratio distortion, TRD), defined as a departure from Mendelian expectations (Crow, 1999), is associated with a wide variety of mechanisms that affect the sperm and ova fertility and viability at developmental stages of the reproductive cycle (e.g., embryos, postnatal metabolism and growth, etc.). Although the genetic background of TRD is mainly related to single-locus factor (i.e., within a single gene or blocks of physically linked loci), interactions between different loci (i.e., epistasis) could also potentially cause specific TRD patterns (Corbett-Detig et al., 2013; Leppälä et al., 2013). Within a TRD context, compatible allele (or genotype) combinations at different loci increase fitness (Fickel and Weyrich, 2011), whereas unfavorable allelic (or genotype) combinations are under-represented, exhibiting deviation from Mendelian proportions (Montagutelli et al., 1996). In addition, as epistasis influences the expressivity of the involved loci, it has been hypothesized that the incomplete penetrance of lethal mutations may involve epistatic interaction (Ballinger and Noor, 2018).

There are documented cases in model organisms, plants and fish where epistatic interactions explain the variation of the observed TRD (Montagutelli et al., 1996; van Boven and Weissing, 2001; Borowsky and Cohen, 2013; Behrouzi and Wit, 2018). Epistatic TRD at different loci has been related to different sources, such as the lethality of embryos carrying particular combinations of alleles (Montagutelli et al., 1996) and lower fertility or diseases (Behrouzi and Wit, 2018), among others. Particularly, epistatic TRD has been related to genic incompatibilities (or hybrid incompatibility) in crosses where deleterious epistatic interactions between heterospecific alleles lead to hybrid sterility and/or inviability (Turelli and Orr, 2000; Fishman et al., 2008; Cutter, 2012; Corbett- Detig et al., 2013; Giesbers et al., 2019; Arends et al., 2022). These genic incompatibilities between diverged genomes, termed Bateson-Dobzhansky-Muller (BDM) incompatibilities, were first described by Bateson (1909), Dobzhansky (1937) and Muller (1942), representing a classic model for the evolution of reproductive isolation in diverging lineages.

Despite the additive paradigm being widely predominant for the technical statistical and computational limitations, the importance of interaction components of the genetic architecture of traits has recently been emphasized (Carlborg and Haley, 2004; Phillips, 2008). It must be highlighted that the relative importance of gene-by-gene interaction has been around for more than a century (Gilbert-Diamond and Moore, 2011) and was recognized as an explanation of the observed deviations from Mendelian ratios (e.g., Bateson, 1909; Neel and Schull, 1954; Hollander, 1955; Phillips, 1998). Identification of epistatic interaction in the genome will provide an opportunity to understand the epistatic effect on agronomically important complex traits allowing insights into the genetic background and complexity underlying reproduction. The TRD approach does not require phenotypic records, but only genomic information in trios (sire-dam-offspring), which is available nowadays for most of the livestock species what opens alternative genomic strategies with potential outcomes to improve reproductive success.

The conventional method of identifying epistatic TRD relies on performing a screen for pairs of loci showing deviations from the expected assuming random segregation and the majority of methods are restricted to an F2 design. The allelic and genotypic TRD models developed by Casellas et al. (2012, 2014, 2017) for single-locus, which are based on tracing allele inheritance from parents to offspring, are flexible and accommodate all types of matings and pedigree structure. Here we adapted these models to account for epistatic interactions and to be applicable for livestock populations. The main objectives of this research, using parent-offspring genotyped trios (sire-dam-offspring) of Holstein cattle, were to investigate the genetic basis for incomplete deviation of single-locus TRD and unraveling pairs of unlinked genomic regions across the whole genome that are not transmitted according to Mendelian inheritance rules, but display epistatic transmission ratio distortion patterns.

2 Materials and methods

2.1 Genotypes and trios

The dataset used to investigate epistatic TRD consisted of 436,651 Holstein genotypes provided by the Canadian Dairy Network database (Lactanet, Guelph, Ontario, Canada). The number of genotyped trios (sire-dam-offspring) was 283,817 and were sampled from all available Holstein genotypes (>1 million; October 2017) with offspring genotyped within 90 days of the birth, thus minimizing selection artifacts (Id-Lahoucine and Casellas, 2017; Id-Lahoucine et al., 2019a). The number of genotyped sire and dams was 5,224 and 117,316, respectively. Animals were genotyped with different SNP genotyping arrays ranging from 6,909 to 777,962 SNPs (Supplementary Material S1) and mapped to the UMD3.1 Bos taurus genome assembly. This data was previously imputed using FImpute (Sargolzaei et al., 2014) to 47,910 SNPs (Id-Lahoucine, 2020).

2.2 Analytical models statistical analyses

Epistatic transmission ratio distortion was evaluated on pairs of SNPs across the whole genome (47,910 SNPs) using three different methods. In order to avoid confounding due to linkage disequilibrium (LD) within chromosomes and to reduce the computational time, epistatic TRD scan was restricted to inter-chromosomal pairs.

Genotypic parameterization of epistatic transmission ratio distortion. Following Casellas et al. (2012, 2020) genotypic TRD model of one single-locus and Kempthorne (1969) epistatic theory, epistatic TRD parameters were included in the probability of the offspring (P_off) for each combination of genotypes of a given pair of two loci (A and B) as:

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{A}\mathrm{B}\mathrm{B}\right)=1+1{\mathrm{\alpha }}_{\mathrm{A}}+1{\mathrm{\alpha }}_{\mathrm{B}}-1{\mathrm{\delta }}_{\mathrm{A}}-1{\mathrm{\delta }}_{\mathrm{B}}+1\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}+1\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}+1\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{a}\mathrm{B}\mathrm{B}\right)=2+0{\mathrm{\alpha }}_{\mathrm{A}}+1{\mathrm{\alpha }}_{\mathrm{B}}+1{\mathrm{\delta }}_{\mathrm{A}}-1{\mathrm{\delta }}_{\mathrm{B}}+0\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}+0\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}+1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{a}\mathrm{a}\mathrm{B}\mathrm{B}\right)=1-1{\mathrm{\alpha }}_{\mathrm{A}}+1{\mathrm{\alpha }}_{\mathrm{B}}-1{\mathrm{\delta }}_{\mathrm{A}}-1{\mathrm{\delta }}_{\mathrm{B}}-1\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}+1\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{A}\mathrm{B}\mathrm{b}\right)=2+1{\mathrm{\alpha }}_{\mathrm{A}}+0{\mathrm{\alpha }}_{\mathrm{B}}-1{\mathrm{\delta }}_{\mathrm{A}}+1{\mathrm{\delta }}_{\mathrm{B}}+0\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}+0\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}+1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{a}\mathrm{B}\mathrm{b}\right)=4+0{\mathrm{\alpha }}_{\mathrm{A}}+0{\mathrm{\alpha }}_{\mathrm{B}}+1{\mathrm{\delta }}_{\mathrm{A}}+1{\mathrm{\delta }}_{\mathrm{B}}+0\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}+0\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}+0\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{a}\mathrm{a}\mathrm{B}\mathrm{b}\right)=2-1{\mathrm{\alpha }}_{\mathrm{A}}+0{\mathrm{\alpha }}_{\mathrm{B}}-1{\mathrm{\delta }}_{\mathrm{A}}+1{\mathrm{\delta }}_{\mathrm{B}}+0\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}+1\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}+0\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}+1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{A}\mathrm{b}\mathrm{b}\right)=1+1{\mathrm{\alpha }}_{\mathrm{A}}-1{\mathrm{\alpha }}_{\mathrm{B}}-1{\mathrm{\delta }}_{\mathrm{A}}-1{\mathrm{\delta }}_{\mathrm{B}}-1\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}+1\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{a}\mathrm{b}\mathrm{b}\right)=2+0{\mathrm{\alpha }}_{\mathrm{A}}-1{\mathrm{\alpha }}_{\mathrm{B}}+1{\mathrm{\delta }}_{\mathrm{A}}-1{\mathrm{\delta }}_{\mathrm{B}}+0\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}+0\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}+1\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}+1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{a}\mathrm{a}\mathrm{b}\mathrm{b}\right)=1-1{\mathrm{\alpha }}_{\mathrm{A}}-1{\mathrm{\alpha }}_{\mathrm{B}}-1{\mathrm{\delta }}_{\mathrm{A}}-1{\mathrm{\delta }}_{\mathrm{B}}+1\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}-1\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}$

where; α_A and δ_A are additive- and dominance-TRD parameters for locus A, respectively, α_B and δ_B are additive- and dominance-TRD parameters for locus B, respectively, αα_e, αδ_e, δα_e and δδ_e are additive-by-additive, additive-by-dominance, dominance-by-additive and dominance-by-dominance epistatic TRD parameters, respectively. The parametric space considered for all the parameters ranges from −1 to 1. Due to the complexity of this model and the number of parameters involved, negative probabilities must be rescaled to 0 for simplification. An additional restriction is needed to guarantee P_off (AABB) + P_off (AaBB) + P_off (aaBB) + P_off (AABb) + P_off (AaBb) + P_off (aaBb) + P_off (AAbb) + P_off (Aabb) + P_off (aabb) = 1. On the other hand, it must be emphasized that all the equations of (offspring) probabilities described here correspond to the unique case of heterozygous-by-heterozygous mating for both loci (i.e., AaBb × AaBb mating), and the probabilities needs to be adapted for each specific mating accordingly. For instance, the P_off (AABB) becomes 4 + 1 α_A + 1 α_B - 1 δ_A - 1 δ_B + 1 αα_e + 1 αδ_e + 1 δα_e - 1 δδ_e whereas P_off (AABb) is 0 for AaBB × AaBB mating. The total number of possible informative matings is 65 and which could be combined in 27 types of matings summarized in the Supplementary Material S2.

Under a Bayesian implementation, the conditional posterior probabilities of the TRD parameters are defined as:

$\mathrm{p}\left({\mathrm{\alpha }}_{\mathrm{A}},{\mathrm{\delta }}_{\mathrm{A}},{\mathrm{\alpha }}_{\mathrm{B}},{\mathrm{\delta }}_{\mathrm{B}},\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}},\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}},\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}},\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}\left|\mathbf{y})\propto \mathrm{p}(\mathbf{y}\right|{\mathrm{\alpha }}_{\mathrm{A}},{\mathrm{\delta }}_{\mathrm{A}},{\mathrm{\alpha }}_{\mathrm{B}},{\mathrm{\delta }}_{\mathrm{B}},\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}},\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}},\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}},\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}\right)\mathrm{p}\left({\mathrm{\alpha }}_{\mathrm{A}}\right)\mathrm{p}\left({\mathrm{\delta }}_{\mathrm{A}}\right)\mathrm{p}\left({\mathrm{\alpha }}_{\mathrm{B}}\right)\mathrm{p}\left({\mathrm{\delta }}_{\mathrm{B}}\right)\mathrm{p}\left(\mathrm{\alpha }{\mathrm{\alpha }}_{\mathrm{e}}\right)\mathrm{p}\left(\mathrm{\alpha }{\mathrm{\delta }}_{\mathrm{e}}\right)\mathrm{p}\left(\mathrm{\delta }{\mathrm{\alpha }}_{\mathrm{e}}\right)\mathrm{p}\left(\mathrm{\delta }{\mathrm{\delta }}_{\mathrm{e}}\right)$

where: y is the column vector of genotypes of the offspring generation.

In order to reduce the computational time, only pairs of SNPs with ≥1,000 informative offspring and ≥20 heterozygous sires or/and ≥50 heterozygous dams were used. Secondly, the observed and expected number of offspring for the 27 combined matings were determined for each pair of loci. The analyses included only pairs of SNPs with ≥1,000 under- or over-represented offspring for a specific two-locus genotype offspring and mating and this difference (between the observed and expected) being relatively 10% higher/lower than the expected number of offspring. Notice that these criteria and cut-off values were chosen based in the results observed in single-locus TRD analyses (Id-Lahoucine, 2020) and all these steps guaranteed a relevant deviation from Mendelian proportions in the corresponding regions allowing to target regions with potential epistatic TRD signals and reducing the computational time considerably. The model was implemented in a Bayesian approach with the metropolis-Hastings (Hastings, 1970) sampling technique. A preliminary scan was implemented with a unique Monte Carlo Markov chain of 11,000 iterations (the first 1,000 were discarded for TRD estimation). The statistical significance of TRD parameters was tested by a Bayes factor (BF, Kass and Raftery, 1995). More accurate estimates were obtained using 550,000 iterations for the selected relevant epistatic regions from the preliminary analyses. In addition, in order to focus mainly on epistasis phenomenon, the ratio of maximum (BF_αα, BF_αδ, BF_δα, BF_δδ)/maximum (BF_αA, BF_δA, BF_αB, BF_δB) > 1,000 was used to select the regions with more statistical evidence of epistatic effects rather direct effects. The BF was also used for the multiple test correction where the top 0.1% significant regions were selected as the most relevant epistatic TRD regions.

Allelic parameterization of epistatic transmission ratio distortion. An allelic epistatic TRD could be targeted by tracing back the inheritance of a combination of two specific alleles of the two loci (i.e., an artificial haplotype) from parents to offspring. Thus, the probability of an offspring (P_off) for each combination of genotypes of a given par of two loci (A and B) was parameterized to include both direct and epistatic allelic TRD effects as:

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{A}\mathrm{B}\mathrm{B}\right)=\left[\left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{A}\mathrm{b}\mathrm{b}\right)=\left[\left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{a}\mathrm{a}\mathrm{B}\mathrm{B}\right)=\left[\left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{a}\mathrm{a}\mathrm{b}\mathrm{b}\right)=\left[\left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{A}\mathrm{B}\mathrm{b},\mathrm{A}\mathrm{A}\mathrm{b}\mathrm{B}\right)=\left[\left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{a}\mathrm{a}\mathrm{B}\mathrm{b},\mathrm{a}\mathrm{a}\mathrm{b}\mathrm{B}\right)=\left[\left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{a}\mathrm{B}\mathrm{B},\mathrm{a}\mathrm{A}\mathrm{B}\mathrm{B}\right)=\left[\left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{a}\mathrm{b}\mathrm{b},\mathrm{a}\mathrm{A}\mathrm{b}\mathrm{b}\right)=\left[\left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{a}\mathrm{B}\mathrm{b},\mathrm{a}\mathrm{A}\mathrm{b}\mathrm{B}\right)=\left[\left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)$

${\mathrm{P}}_{\mathrm{o}\mathrm{f}\mathrm{f}}\left(\mathrm{A}\mathrm{a}\mathrm{b}\mathrm{B},\mathrm{a}\mathrm{A}\mathrm{B}\mathrm{b}\right)=\left[\left(0.5+{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5+{\mathrm{\beta }}_{\mathrm{B}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{A}}\right)\times \left(0.5-{\mathrm{\beta }}_{\mathrm{B}}\right)\right]\times \left(1+{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)\times \left(1-{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)$

where; β_A and β_B are direct TRD parameters for locus A and B, respectively, β_AB/Ab, β_AB/aB, β_AB/ab, β_Ab/aB, β_Ab/ab and β_aB/ab are the 6 heterozygous pairwise combinations for epistatic TRD parameters. Note that AB, Ab, aB and ab are the 4 possible artificial haplotype alleles of the two implicated SNPs. For all parameters, flat priors were assumed within a parametric space ranging from −0.5 to 0.5 for direct TRD effects and from −1.0 to 1.0 for epistatic TRD effects. The parametric space of direct effect is based on the principles of Mendelian inheritance (i.e., the probability of transmission of one specific allele ranges from 0 (e.g., β_A = −0.5) to 1 (e.g., β_A = 0.5), and where 0.5 (e.g., β_A = 0.0) corresponds to null TRD). In contrast, the range from −1 to 1 for epistatic parameters could be viewed as the positive/negative preferential transmission of one combination of two alleles (artificial haplotype) or its opposite heterozygous combinations, respectively, whereas 0 indicates null epistatic TRD. Also, these probabilities corresponded to the unique case of heterozygous-by-heterozygous mating (i.e., AaBb × AaBb), and the probabilities must be adapted for each specific mating accordingly.

Under a Bayesian implementation, the conditional posterior probabilities of the TRD parameters are defined as:

$\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{A}},{\mathrm{\beta }}_{\mathrm{B}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{A}\mathrm{b}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{B}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{b}},{\mathrm{\beta }}_{\mathrm{a}\mathrm{B}/\mathrm{a}\mathrm{B}}\left|\mathbf{y})\propto \mathrm{p}(\mathbf{y}\right|{\mathrm{\beta }}_{\mathrm{A}},{\mathrm{\beta }}_{\mathrm{B}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{A}\mathrm{b}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{B}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}},{\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{b}},{\mathrm{\beta }}_{\mathrm{a}\mathrm{B}/\mathrm{a}\mathrm{B}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{A}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{B}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{A}\mathrm{b}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{B}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{A}\mathrm{B}/\mathrm{a}\mathrm{b}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{B}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{A}\mathrm{b}/\mathrm{a}\mathrm{b}}\right)\mathrm{p}\left({\mathrm{\beta }}_{\mathrm{a}\mathrm{B}/\mathrm{a}\mathrm{B}}\right)$

where; y is the column vector of genotypes of the offspring generation.

This method is similar to the heterozygous pairwise combination procedure of haplotype analysis (single locus) described by Id-Lahoucine et al. (2019a), which is highly computationally demanding. For this reason, a simplified algorithm was implemented for a preliminary scan. The biallelic-haplotype procedure for haplotype analysis described by Id-Lahoucine et al. (2019a) was adapted to estimate TRD for the artificial haplotypes in order to target signals of epistasis. Thus, the transmission probability (P) of an artificial haplotype from heterozygote parents to offspring was parameterized including one overall epistatic TRD effect (α_j) for each specific j artificial haplotype (AH) allele (P (AH_j) = 1—P (AH_-j) = 0.5 + α_j).

As an exhaustive search for epistatic allelic TRD is computationally unfeasible, epistatic overall TRD were estimated only for the combination of pairs of SNPs which include one SNP with significant direct TRD effect already identified from analysis of single-locus TRD (2,962 SNPs; Id-Lahoucine, 2020). After the construction of artificial bi-allelic haplotypes, the analyses were performed within a Bayesian framework using a TRDscan v.1.0 software (Id-Lahoucine et al., 2019a) with a unique Monte Carlo Markov chain of 110,000 iterations where the first 10,000 iterations were discarded as burn-in. The statistical significance of TRD was evaluated using a Bayes factor (Kass and Raftery, 1995). Following this analysis, the candidate regions with signals of epistatic allelic pattern obtained by the simplified method were re-examined with the full epistatic allelic model (with 550,000 iterations), thus, allowing to overcome the limitation of this simplified method which was not parameterized to take into account both direct and epistatic effects simultaneously. Finally, both allelic and genotypic parameterizations were compared using the deviance information criterion (DIC, Spiegelhalter et al., 2002) to determine the goodness-of-fit and the epistatic pattern of each pair of regions.

Double recessive epistatic transmission ratio distortion. The criteria applied for both genotypic and allelic model will discard the possibility of detecting double recessive epistatic TRD. In this sense, in order to target double recessive epistatic TRD, the biallelic-haplotype procedure with artificial haplotypes was implemented to screen for absent of double homozygous offspring for two-locus genotypes. For this, a minimal 15 none-observed offspring for one single combination of homozygous genotypes were chosen, while ensuring no deviation from Mendelian proportions for other genotypes.

Multiple test correction: For multiple test correction, the BF was used in order to choose the best candidate regions with epistatic effects. Thus, the multiple test correction was performed by selecting the top 0.1% of the regions within each category according to BF.

3 Results

3.1 Genotypic epistatic transmission ratio distortion

Pairs of regions with epistatic TRD estimates were found widely across the Holstein genome. Using a threshold of ≥1,000 for BF, a large number of pairs of SNPs were observed with decisive evidence for epistatic TRD. In order to discard false TRD and to identify the most relevant regions, the following steps were considered. Firstly, in order to focus on epistasis phenomenon, regions with more statistical relevance for direct effects than epistatic effects were discarded. Thus, only region with a minimal ratio of 1,000 between the maximum BF of epistatic and direct effects were considered (i.e., maximum (BF_αα, BF_αδ, BF_δα, BF_δδ)/maximum (BF_αA, BF_δA, BF_αB, BF_δB) > 1,000). Secondly, estimates with large credible intervals were discarded as artifacts of the convergence of the model following Id-Lahoucine (2020). For this purpose, pairs of SNPs with a coefficient of variation >20% for significant epistatic effects were excluded. Notice that from single-locus TRD, when clear TRD exist, the standard variation of TRD estimates is mostly null (<0.01), thus, the use of coefficient of variation of TRD estimates is a straightforward rule to discard regions with instable convergence. Thirdly, given that most regions have several direct and epistatic effects with different statistical significance, the results were separated into their corresponding effects (i.e., additive-by-additive, additive-by-dominance/dominance-by-additive or dominance-by-dominance) according to the most relevant effect in terms of BF. After filtering the results following the previous steps, the number of the obtained pairs were 59,831, 87,699 and 20,549 for additive-by-additive, additive-by-dominance/dominance-by-additive and dominance-by-dominance effects, respectively.

After implementing a multiple test correction, which was based on selecting the top 0.1% of the regions within each category according to BF, the number of regions reduced to 169, being 60, 88 and 21 with additive-by-additive, additive-by-dominance/dominance-by-additive and dominance-by-dominance effects, respectively. This strategy was implemented with the objective to select the most significant “top hits” of epistatic TRD signals following François et al. (2016) and Boschiero et al. (2018). These pairs of regions were chosen for a second analysis with a large Monte Carlo Markov chain of 500,000 iterations. Applying the previous criteria described above to the new estimated, 55 out of 88 regions with additive-by-dominance/dominance-by-additive were discarded. Moreover, among the identified pairs, we observed that some SNPs interacted with more than one individual SNP. For instance, the SNP BTA17:39,696,262 interacted with 26 SNPs covering from 4,197,354 bp to 12,447,484 bp on BTA23. Within this context, only the pair of SNPs with the highest BF were maintained as the best candidate regions when one single locus interacts with several physically linked SNPs. Thus, the number of regions reduced to 7, 19 and 6 with additive-by-additive, additive-by-dominance/dominance-by-additive and dominance-by-dominance effects, respectively. These results are summarized in Table 1 for additive-by-additive and dominance-by-dominance effects and in Table 2 for additive-by-dominance/dominance-by-additive. Additional details (e.g., number of informative parents, frequencies, etc.) are presented in the Supplementary Material S3.

TABLE 1

TABLE 1. Pairs of genomic regions identified with epistatic transmission ratio distortion (TRD) in Holstein cattle with mainly additive-by-additive or dominance-by-dominance effects.

TABLE 2

TABLE 2. Pairs of genomic regions identified with epistatic transmission ratio distortion (TRD) in Holstein cattle with mainly additive-by-dominance effects.

When comparing between the genotypic and allelic models, different goodness-of-fit values were observed among the regions with additive-by-additive effects displaying reductions from 1,186.35 to 46,445.15 DIC units favoring the genotypic model over the allelic model. Notice that models with smaller DIC values indicate a better fit, and differences between models greater than 3 DIC units are considered statistically relevant (Spiegelhalter et al., 2002). In the case of regions with dominance-by-dominance effects, reductions up to 49,728.74 DIC units were observed. In addition, out of the regions found with additive-by-dominance/dominance-by-additive, only 8 pairs of regions favored the genotypic model with differences of up 20,349.96 DIC units. Therefore, the remaining pairs of regions (11) detected with the genotypic model with additive-by-dominance/dominance-by-additive showed better goodness-of-fit with the allelic model in terms of DIC. Nevertheless, more statistical relevance based on BF was observed for TRD parameters of genotypic model (up to log₁₀(BF) = 1,964.13) in comparison to TRD estimated of allelic model (up to log₁₀(BF) = 269.55). It is important to remember that DIC was computed based on all TRD parameters combined and the regions here were separated and selected based on their type of effects (e.g., additive-by-additive, additive-by-dominance) and BF.

3.2 Allelic epistatic transmission ratio distortion

Using the simplified method of the allelic model, a minimal increase of 0.05 in the TRD magnitude explained by the artificial haplotype (i.e., gamete) compared to the single SNP and an equal or greater number of under-/over-represented offspring (i.e., ≈ α_j × 2 × number of informative offspring) were considered to identify candidate regions with allelic epistatic TRD pattern. The total number of obtained regions detected with BF ≥ 1,000 was 4,852. For the identified artificial haplotypes, the absolute TRD magnitude ranged from |0.06| to |0.42|, and the maximum number of under- and/or over-represented offspring was 12,756 genotypes. The additional magnitude of TRD explained by artificial haplotypes reached a maximum of 0.25, with 8, 166 and 569 pairs with 0.20, 0.15 and 0.10, respectively. For example, pairs from the top regions (with highest BF) were found with individual SNP (BTA23:6,948,746) displaying an overall TRD of −0.22 and 3,078 under-represented offspring (log₁₀(BF) = 302.62) and exhibited a deviation of −0.37 and 3,908 under-represented offspring (log₁₀(BF) = 692.38) when paired with a specific allele in another SNP (BTA1:45,157,959). Thus, when the same SNP was found interacting in several pairs, we considered the pair of loci with the highest BF as the best candidate region with epistatic TRD. This reduced the number of pairs of regions to 78 (Supplementary Material S4). In terms of goodness-of-fit, among the 78 pairs of regions detected with the allelic model, minimal DIC values were observed on 76 pairs when compared to the genotypic model. Specifically, reduction from 19.38 up to 23,647.43 DIC units were observed for the allelic model in comparison to the genotypic model. The re-analyses of these regions with the full allelic model, with epistatic and direct TRD effects, displayed significant epistatic TRD effect with log₁₀(BF) > 9.61, being log₁₀(BF) = 969.06 the maximum value observed for an epistatic heterozygous pairwise combination. The number of regions that showed clear epistatic effects with the full model in all heterozygous pairwise combinations including the artificial haplotype were 29 and reduced to 17 when only one single artificial haplotype displayed epistatic TRD (Table 3).

TABLE 3

TABLE 3. Pairs of genomic regions identified with epistatic transmission ratio distortion (TRD) in Holstein cattle with allelic patterns.

3.3 Recessive epistatic transmission ratio distortion

The total number of pairs of regions with at least 15 expected homozygous offspring but none of them observed assuming random assortment for both implicated SNPs, was 67. This was reduced to 56 after discarding pairs of regions pointing to physically related regions (Table 4, Supplementary Material S5). The maximum number of expected homozygous offspring for both SNPs, but none of them observed, was 23. The number of informative sires, dams and offspring for the identified regions reached to 105, 2,809 and 11,995, respectively.

TABLE 4

TABLE 4. Pairs of genomic regions identified with epistatic transmission ratio distortion (TRD) in Holstein cattle with recessive patterns.

4 Discussion

The current accessibility of high-throughput genotyping technologies in big data era have allowed the investigation of epistasis more deeply. Previous studies used genotype data and phenotypes or gene expression to investigate epistatic interactions between loci (e.g., Strange et al., 2013; Hemani et al., 2014; Huang et al., 2014; Mackay, 2014). Nevertheless, though many efforts have been done for identifying epistatic interactions, the methods/algorithmics are still being developed (e.g., Beissinger et al., 2015; Sun et al., 2016; 2017; Behrouzi and Wit, 2018; Id-Lahoucine et al., 2019b). In this paper we extended TRD models from single-locus to capture epistatic interactions, presenting a practical methodology when trios of parent-offspring genotypes (sire-dam-offspring) are available, not being restricted to only F2 designs, which is the most commonly used design in recent studies of TRD and epistasis in model organisms (e.g., Brennan et al., 2014; Niedzicka et al., 2017; Haddad et al., 2018).

4.1 Genotypic epistatic transmission ratio distortion

The Bayes factor is the standard Bayesian tool to compare two competing models (Kass and Raftery, 1995). For determining statistical significance, the possible choices of a BF threshold could be variable, BF ≥ 10 indicating strong evidence and BF ≥ 100 more decisive evidence according to the Jeffreys’ (1984) scale. Within the context of TRD analysis, the BF can be viewed as a measure of the strength of statistical evidence of the detected regions, given both TRD magnitude and sample size of informative offspring simultaneously. The value of BF increases with the increase of TRD magnitude and the number of informative offspring (Id-Lahoucine et al., 2019a). Here, a higher threshold of BF (≥1,000) was considered to ensure targeting epistatic regions with relevant TRD magnitudes. Even though the prevalence of signals of epistatic TRD widely extended across the Holstein genome (i.e., a large number of pairs of SNPs (>150,000) were initially identified), it is expected that part of them are simply artifacts of the sampling fluctuations that generate random/false TRD. It should be noted that the complexity and dimensionality of the model itself with 8 parameters and their possible interactions in the same data could generate false epistatic TRD. It was previously reported that the structure of the data with an unbalanced number of trios across matings could easily generate TRD artifacts (Casellas et al., 2020; Id-Lahoucine, 2020). For example, an inappropriate convergence of TRD estimation is produced when combinations of different TRD parameters maximize the likelihood of the data and result in large credible intervals for TRD estimates (including often zero values). For this reason, a minimal dispersion of TRD effects was assumed to discard potential TRD artifacts.

Given the scope of the study, which focus mainly on epistasis phenomenon, the ratio of maximum (BF_αα, BF_αδ, BF_δα, BF_δδ)/maximum (BF_αA, BF_δA, BF_αB, BF_δB) was used to select the regions with more statistical evidence of epistatic effects rather direct effects. Assuming, for example, that BF_αα [= p (αα_e = 0)/p (αα_e = 0|y)] and BF_αA [= p (α_A = 0)/p (α_A = 0|y)] are the maximum values, this ratio becomes = p (α_A = 0|y)/p (αα_e = 0|y). Given p (α_A = 0|y) = (p (y|α_A = 0)p (α_A = 0))/p(y), the ratio simplifies to p (y|α_A = 0)/p (y|αα_e = 0), a simple and direct measure of comparison between direct and epistatic effects. Thus, by using a ratio of 1,000, we ensured that the epistatic effect is much more relevant to explain the observed data of the implicated regions than the direct effects.

After the multiple test correction (top 0.1% of the regions within each category according to BF), the minimal and maximal log₁₀(BF) obtained after the multiple test correction was 1,223.8 and 20,223, respectively. Notice that BF is a measure of the strength of statistical evidence (for both TRD magnitude and sample size of informative offspring). The likelihood of the data assuming epistatic effects for these regions is ≥ 10^1223.8 times more probable that assuming null TRD effects, very strongly supporting the presence of epistatic effects.

It is important to mention that while the number of regions with additive-by-additive and dominance-by-dominance did not reduce after using large sample chains, they were considerably reduced for regions with additive-by-dominance/dominance-by-additive effects. This latter observation is partially due to the confounding between the direct dominance effect and the additive-by-dominance/dominance-by-additive epistatic effects as observed in different estimates among different chains in some cases. This suggests that several combinations of TRD estimates maximize the likelihood of the data, as mentioned before, and therefore, it is more difficulty to estimate their effects. Moreover, regions with additive-by-additive effects showed the clearest patterns of epistatic effects. One of the regions with an estimated αα_e = 0.99 (standard deviation = 0.00015, log₁₀(BF_αα) = 4,880.73) displayed an over-representation of 8,054 offspring (40.04%; 20,117 observed vs. 12,063 expected) for AABB and 14,800 (30.29%; 48,867 observed vs. 34,067 expected) for aaBB, whereas an under-representation of 5,567 (49.80%; 5,633 observed vs. 11,179 expected) for aaBB and 4,227 (66.85%; 2,100 observed vs. 6,323 expected) for AAbb (Figure 1A, Supplementary Material S3). On the other hand, regarding the dominance-by-dominance effect, a pair of regions with δδ_e = −0.99 (standard deviation = 0.00004, log₁₀(BF_δδ) = 16,219.5) showed an under-representation of 14,338, 14,079, 25,459 and 24,814 for AaBB, AABb, aaBb and Aabb, respectively, and an over-representation of 17,864, 1, 37,346, 9 and 29,853 for AABB, aaBB, AaBb, AAbb and aabb, respectively, following the parameterization for dominance-by-dominance effects (Figure 1B, Supplementary Material S3). Those examples are the simplest cases with clear epistatic TRD patterns, but it must be taken into consideration that most regions have several direct and epistatic effects with different statistical significance as evidence of the complexity of the epistatic phenomenon in the genome of cattle.

FIGURE 1

FIGURE 1. The observed, under- and over-represented offspring across matings in two-locus genotype with (A) additive-by-additive and (B) dominance-by-dominance epistatic effects.

Furthermore, identification of SNPs interacting with several physically linked SNPs supports the relevance of the epistasis phenomenon. This observation is potentially explained by the LD, where one individual or different physically linked SNPs (i.e., a genomic region/locus) interact with one or different physically linked SNPs, supporting the importance of the epistatic interaction found. For this reason, only the pair of SNPs with the highest BF were maintained as the best candidate regions when one single locus interacted with several physically linked SNPs. On the other hand, an individual SNP interacting with several not linked SNP were also found. This latter case might potentially suggest a multiple interaction involving different unlinked loci (of order 3 or more).

4.2 Allelic epistatic transmission ratio distortion

The preliminary scan with the allelic model was implemented only in pairs of loci with already identified SNP with direct TRD effect in order to target possible epistatic interaction explaining the incomplete penetrance (deviation) of single-locus TRD. The criteria implemented in the simplified method were used to ensure more TRD was explained with epistatic interaction and that the observed TRD of the individual SNP was captured by an artificial haplotype. Explicitly, the restriction of an equal or greater number of under-/over-represented offspring was used to guarantee that the under-/over-represented offspring generated from a specific parent with a particular SNP allele corresponds to a single artificial haplotype (gamete) on those parents. Note that the linkage of specific alleles in a single gamete could be a result of selection or random effects, as it is known by theory that selection or genetic drift could favor combinations of alleles and consequently induce LD at physically unlinked loci (Bulmer, 1971; Ohta, 1982). However, by using trios of parent-offspring genotypes with offspring genotyped at early age, the issue of selection potentially was minimized, and the detected regions are expected to be associated with the reproduction cycle. In this case, under the hypothesis that the partial lethality of one allele at the first locus is linked to a specific allele on the second paired locus, this methodology could allow more insights to be gained in understating the penetrance of lethal alleles when knowing the genotype of the implicated loci. For instance, the variation of penetrance explained by the interaction could be due to some functional interaction of some unfavorable allelic combinations when inherited together in the genome of the progeny and probably affecting the fertility of parents and/or viability of embryos/offspring.

On the other hand, it is important to mention that, as the analyses are restricted to SNPs with single-locus TRD in which the majority had low frequencies (Id-Lahoucine, 2020), the multiple test combination with all other SNPs could easily generate allelic epistasis in some of them by chance. Indeed, SNPs interacting with several SNPs were observed. Specifically, only 11 pairs of SNPs had one single interaction whereas the maximum number of SNPs observed interacting with one individual SNP was 866. Even though it is plausible to assume that there are multiple interactions between several loci, these multiple interactions must be taken with caution given the low frequencies of the analyzed SNPs (from the 2,962 used SNPs, 1,424 and 950 had frequency <0.05 and 0.01, respectively). However, given the analysis were restricted to candidate regions (with low frequencies in most cases), the pair of loci with the highest BF was considered as the best candidate region with epistatic TRD, when the same SNP was found interacting with several SNPs.

The results of the analysis using the full model including simultaneously both direct and epistatic effects supported the previous findings. First, this emphasizes the effectiveness of the simplified approach to largely reduce the computational time. Second, in comparison to the simplified method, the epistatic effects observed in one artificial haplotype were also observed in all (3) or some of the heterozygous pairwise combinations including the artificial haplotype identified. An example of this was an artificial haplotype (aB) with a TRD magnitude of 0.37 (log₁₀(BF) = 769.12) that, using biallelic-haplotype method, showed a TRD of −0.78 (log₁₀(BF) = 351.96), −0.67 (log₁₀(BF) = 216.05) and 0.74 (log₁₀(BF) = 192.72) for β_AB/aB, β_Ab/aB and β_aB/ab, respectively, whereas keeping β_A = −0.31 (log₁₀(BF) = 471.30) and β_B = 0.29 (log₁₀(BF) = 328.74) for direct effect. The number of regions which showed clear epistatic effects in all heterozygous pairwise combinations including the artificial haplotype, was 29. Nevertheless, in some cases, epistatic TRD was also observed in heterozygous pairwise combinations with other artificial haplotype alleles. This result adds difficulties in understanding epistatic TRD, but could be explained by the different linkage with the causal mutations among families, higher order of epistasis interactions, and does not rule out the possibility of involvement of specific interactions between the artificial haplotypes. A total of 17 regions were identified when only pairs of regions displaying epistatic TRD for a single artificial haplotype were considered (Table 3).

Moreover, we observed that the magnitude obtained from single-locus analyses were similar (mostly identical) to those obtained from the epistatic allelic model (direct effects) when both direct and epistatic TRD effect were simultaneously estimated. This interesting behavior support the robustness and the accuracy of allelic model, where basically the model it tried to track the inheritance of alleles from parents to offspring independently of its source. It is also important to mention that a similar behavior was already reported in single-locus TRD analyses, thus highlighting the robustness of allelic model when comparing both allelic and genotypic model (Id-Lahoucine, 2020; Id-Lahoucine et al., 2022). Thus, this latter suggests that epistatic TRD could generate signals of single-locus TRD (as expected) and consequently could make difficult to differentiate between the origins of single-locus TRD in some cases. However, it is important to emphasize that the likelihood ratio between the full and null TRD models exhibited values up to 10^1448.11 (minimal observed was 10^12.71). More specifically, the epistatic parameters seems to be very importance as the likelihood ratio between the full model and an alternative model without epistatic parameters exhibited values up to 10^1,261.89 (minimal 2,630.27), supporting the importance of the epistatic TRD phenomenon in the data.

Finally, it is important to highlight that not always a disadvantage was observed when a specific artificial haplotype was transmitted, but also artificial haplotype alleles with positive effects were found. In this sense, one single artificial haplotype with a negative TRD was observed, associated with an under-represented offspring and potentially lethality of the offspring of carrier parent, whereas the alternative artificial haplotypes had positive effects. In contrast, when one single artificial haplotype had a positive effect (whereas negative or null for the remaining artificial haplotypes), it suggests the advantage of one haplotype to be transmitted over other haplotypes. This could be due to an increase in fitness when an allele is paired with a specific allele at another locus. In addition, it could be hypothesized that these regions with positive selection are more related to biological processes acting before fertilization, where gametes with specific allele combinations have a preferential transmission to the next-generation.

4.3 Recessive epistatic transmission ratio distortion

Absence (or depletion) of individuals with two-locus genotypes in the homozygous statue were found across the Holstein genome, suggesting multiple variants with recessive epistatic TRD pattern. The number of candidate pairs of loci were 56, with 23 of them having the maximum number of expected homozygous offspring, but none of them being observed. It is important to mention that Holstein haplotype 3 was identified initially by VanRaden et al. (2011) with only 7 non-observed homozygous offspring from heterozygous sires in combination with heterozygous maternal grandsires. These findings provide evidence of regions carrying deleterious mutations that are expressed when presented together in a homozygous state, potentially producing lethality.

4.4 Quantitative trait loci and epistatic transmission ratio distortion

To support TRD findings, several quantitative trait loci (QTL) from the Cattle QTL database (CattleQTLdb, www.animalgenome.org/cgi-bin/QTLdb/BT/index) were found to overlap with epistatic TRD regions. Among them, the pair of SNPs BTA1:36209316 and BTA21:26405185 were reported with QTL for reproduction traits such as gestation length (Maltecca et al., 2008) and first service to conception traits (Kiser et al., 2019). Cole et al. (2011) reported QTL for both Calving ease and Stillbirth for BTA23:1674622 and BTA28:2233574. Other pairs of SNPs (BTA 2:27724930 and BTA 7:10835967) were also found with QTL for calving ease and stillbirth by Cole et al. (2011); Müller et al. (2017). In addition, QTL for the interval from first to last insemination and non-return rate (Müller et al., 2017), interval to first estrus after calving (Liu et al., 2017) and conception rate at first service (Galliou et al., 2020) were also found for BAT4:85759993 and BAT21:37978844. For the particular pair “BAT7:10835967 and BAT8:101251865” we find the existence of QTL for the interval from first to last insemination (Zhang et al., 2019) and daughter pregnancy rate (Cole et al., 2011). The SNP BAT7:10835967 also interacted with BAT10:29761646 and where QTL were reported for calving ease, daughter pregnancy rate and stillbirth by Cole et al. (2011), calving ease by Müller et al. (2017), interval to first estrus after calving by Liu et al. (2017) and conception rate at first service by Galliou et al. (2020). Moreover, the same SNP BAT7:10835967 interacted with BAT24:33863680 and both presented a QTL for stillbirth by Thomasen et al. (2008) and Müller et al. (2017). All the results supported the reliability of epistatic TRD findings, however, further research investigating more in-depth the genomic regions with epistatic TRD concerning their biological and functional implications are needed.

5 Conclusion

Our results aimed to elucidate the prevalence and patterns of epistatic transmission ratio distortion across the Holstein genome. Different epistatic TRD patterns were observed. Using genotypic models, 7, 19 and 6 pairs of SNPs with additive-by-additive, additive-by-dominance/dominance-by-additive and dominance-by-dominance effects were identified with decisive evidence, respectively. The allelic TRD model revealed 17 pairs of SNPs that offered more insights into understating the penetrance of single-locus lethal alleles, providing a more exact probability of lethality when the genotype of both implicated loci is taken into consideration. Scanning for the depletion of individuals carrying double homozygous genotypes for unlinked loci compared to expected, revealed 56 pairs of SNPs with a recessive epistatic TRD pattern. Additionally, the detection of epistatic TRD on pairs of SNPs with QTL for reproductive traits supported the reliability of epistatic TRD findings. Finally, in this study, we demonstrated candidate genomic regions harboring epistatic interactions with potential biological implications in economically important traits, such as reproduction.

Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: Data that support the findings of this study are available from the Canadian Dairy Network (a member of Lactanet, Guelph, Ontario, Canada) upon reasonable request to the corresponding author, but restrictions apply to the availability of these data, which were used under a license of a material transfer agreement for the current study, and thus are not publicly available. Requests to access these datasets should be directed to Angela Cánovas, YWNhbm92YXNAdW9ndWVscGguY2E=.

Ethics statement

Ethical review and approval was not required for the animal study because the study was exempt from ethical approval procedures.

Author contributions

SI-L, JC, and AC conceived the research and developed the models. SI-L developed the program for data analyses and performed the corresponding analyses. SI-L wrote the first draft of the manuscript. SI-L, JC, AC, FS, and FM participated in the interpretation of the results and revision of the manuscript. AC and JC supervised the research. AC secured the funding for the project. All authors read and approved the final manuscript.

Funding

Authors thank the financial support in main part by Ontario Ministry of Agriculture, Food, and Rural Affairs (OMAFRA; Ontario, Canada), the Ontario Agri-Food Innovation Alliance, Agriculture and Agri-Food Canada, and by the Sustainable Beef and Forage Science Cluster (FDE.13.17) funded by the Canadian Beef Cattle Check-Off, Beef Cattle Research Council (BCRC), Alberta Beef Producers, Alberta Cattle Feeders’ Association, Beef Farmers of Ontario, La Fédération des Productuers de bovins du Québec, and Agriculture and Agri-Food Canada’s Canadian Agricultural Partnership. This study was also supported by NSERC (Natural Sciences and Engineering Research Council).

Acknowledgments

The authors acknowledge the data provided by Lactanet (Guelph, Ontario, Canada).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1132796/full#supplementary-material

References

Arends, D., Karst, S., Heise, S., Korkuc, P., Hesse, D., and Brockmann, G. A. (2022). Transmission distortion and genetic incompatibilities between alleles in a multigenerational mouse advanced intercross line. GENETICS 220 (1), iyab192. doi:10.1093/genetics/iyab192

PubMed Abstract | CrossRef Full Text | Google Scholar

Ballinger, M. A., and Noor, M. A. F. (2018). Are lethal alleles too abundant in humans? 2018. Trends Genet. 34, 87–89.

PubMed Abstract | CrossRef Full Text | Google Scholar

Bateson, W. (1909). “Heredity and variation in modern lights,” in Darwin and modern science. Editor A. C. Seward (Cambridge, MA, United Kingdom: Cambridge University Press), 85–10.

Google Scholar

Behrouzi, P., and Wit, E. C. (2018). Detecting epistatic selection with partially observed genotype data by using copula graphical models. J. R. Stat. Soc. C 68, 141–160. doi:10.1111/rssc.12287

CrossRef Full Text | Google Scholar

Beissinger, T. M., Gholami, M., Erbe, M., Weigend, S., Weigend, A., de Leon, N., et al. (2015). Using the variability of linkage disequilibrium between subpopulations to infer sweeps and epistatic selection in a diverse panel of chickens. Heredity 116, 158–166. doi:10.1038/hdy.2015.81

PubMed Abstract | CrossRef Full Text | Google Scholar

Borowsky, R., and Cohen, D. (2013). Genomic consequences of ecological speciation in Astyanax cavefish. PLoS ONE 8 (11), e79903. doi:10.1371/journal.pone.0079903

PubMed Abstract | CrossRef Full Text | Google Scholar

Boschiero, C., Moreira, G. C. M., Gheyas, A. A., Godoy, T. F., Gasparin, G., Mariani, P. D. S. C., et al. (2018). Genome-wide characterization of genetic variants and putative regions under selection in meat and egg-type chicken lines. BMC Genomics 25 (1), 83. doi:10.1186/s12864-018-4444-0

CrossRef Full Text | Google Scholar

Brennan, A. C., Hiscock, S. J., and Abbott, R. J. (2014). Interspecific crossing and genetic mapping reveal intrinsic genomic incompatibility between two Senecio species that form a hybrid zone on Mount Etna, Sicily. Heredity 113, 195–204. doi:10.1038/hdy.2014.14

PubMed Abstract | CrossRef Full Text | Google Scholar

Bulmer, M. G. (1971). The effect of selection on genetic variability. Am. Nat. 105, 201–211. doi:10.1086/282718

CrossRef Full Text | Google Scholar

Carlborg, O., and Haley, C. S. (2004). Epistasis: Too often neglected in complex trait studies? Nat. Rev. Genet. 5, 618–625. doi:10.1038/nrg1407

PubMed Abstract | CrossRef Full Text | Google Scholar

Casellas, J., Cañas-Álvarez, J. J., González-Rodríguez, A., Puig-Oliveras, A., Fina, M., Piedrafita, J., et al. (2017). Bayesian analysis of parent-specific transmission ratio distortion in seven Spanish beef cattle breeds. Anim. Genet. 48, 93–96. doi:10.1111/age.12509

PubMed Abstract | CrossRef Full Text | Google Scholar

Casellas, J., Gularte, R. J., Farber, C. R., Varona, L., Mehrabian, M., Schadt, E. E., et al. (2012). Genome scans for transmission ratio distortion regions in mice. Genetics 191, 247–259. doi:10.1534/genetics.111.135988

PubMed Abstract | CrossRef Full Text | Google Scholar

Casellas, J., Id-Lahoucine, S., and Cánovas, A. (2020). Discriminating between allele- and genotype-specific transmission ratio distortion. Anim. Genet. 51 (6), 847–854. doi:10.1111/age.13007

PubMed Abstract | CrossRef Full Text | Google Scholar

Casellas, J., Manunza, A., Mercader, A., Quintanilla, R., and Amills, M. (2014). A flexible Bayesian model for testing for transmission ratio distortion. Genetics 198, 1357–1367. doi:10.1534/genetics.114.169607

PubMed Abstract | CrossRef Full Text | Google Scholar

Cole, J. B., Wiggans, G. R., Ma, L., Sonstegard, T. S., LawlorJr, T. J., Crooker, B. A., et al. (2011). Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows. BMC Genomics 12, 408. doi:10.1186/1471-2164-12-408

PubMed Abstract | CrossRef Full Text | Google Scholar

Corbett-Detig, R., Zhou, J., Clark, A. G., Hartl, D. L., and Ayroles, J. F. (2013). Genetic incompatibilities are widespread within species. Nature 504, 135–137. doi:10.1038/nature12678

PubMed Abstract | CrossRef Full Text | Google Scholar

Crow, J. F. (1999). Unmasking a cheating gene. Science 283, 1651–1652. doi:10.1126/science.283.5408.1651

PubMed Abstract | CrossRef Full Text | Google Scholar

Cutter, A. D. (2012). The polymorphic prelude to Bateson-Dobzhansky-Muller incompatibilities. Trends Ecol. Evol. 27, 209–218. doi:10.1016/j.tree.2011.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobzhansky, T. H. (1937). Genetics and the origin of species. New York, NY: Columbia University Press.

Google Scholar

Fickel, J., and Weyrich, A. (2011). “Female mate choice in rodents,” in From genes to animal behavior. Editors M. Inoue-Murayama, S. Kawamura, and A. Weiss (Japan: Springer Japan), 3–33. doi:10.1007/978-4-431-53892-9_1

CrossRef Full Text | Google Scholar

Fishman, L., Aagaard, J., and Tuthill, J. C. (2008). Toward the evolutionary genomics of gametophytic divergence: Patterns of transmission ratio distortion in monkeyflower (mimulus) hybrids reveal a complex genetic basis for conspecific pollen precedence. Evolution 62, 2958–2970. doi:10.1111/j.1558-5646.2008.00475.x

PubMed Abstract | CrossRef Full Text | Google Scholar

François, O., Martins, H., Caye, K., and Schoville, S. D. (2016). Controlling false discoveries in genome scans for selection. Mol. Ecol. 25 (2), 454–469. doi:10.1111/mec.13513

PubMed Abstract | CrossRef Full Text | Google Scholar

Galliou, J. M., Kiser, J. N., Oliver, K. F., Seabury, C. M., Moraes, J. G. N., Burns, G. W., et al. (2020). Identification of loci and pathways associated with heifer conception rate in U.S. Holsteins. U.S. holsteins. Genes. 11, 767. doi:10.3390/genes11070767

PubMed Abstract | CrossRef Full Text | Google Scholar

Giesbers, A. K. J., den Boer, E., Ulen, J. J. W. E. H., van Kaauwen, M. P. W., Visser, R. G. F., Niks, R. E., et al. (2019). Patterns of transmission ratio distortion in interspecific Lettuce hybrids reveal a sex-Independent gametophytic barrier. Genetics 211, 263–276. doi:10.1534/genetics.118.301566

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilbert-Diamond, D., and Moore, J. H. (2011). Analysis of gene-gene interactions. Curr. Protoc. Hum. Genet. 1-14.

CrossRef Full Text | Google Scholar

Haddad, R., Meter, B., and Ross, J. A. (2018). The genetic architecture of intra-species hybrid mito-nuclear epistasis. Front. Genet. 9, 481. doi:10.3389/fgene.2018.00481

PubMed Abstract | CrossRef Full Text | Google Scholar

Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57, 97–109. doi:10.1093/biomet/57.1.97

CrossRef Full Text | Google Scholar

Hemani, G., Shakhbazov, K., Westra, H. J., Esko, T., Henders, A. K., McRae, A. F., et al. (2014). Detection and replication of epistasis influencing transcription in humans. Nature 508, 249–253. doi:10.1038/nature13005

PubMed Abstract | CrossRef Full Text | Google Scholar

Hollander, W. F. (1955). Epistasis and hypostasis. J. Hered. 46, 222–225. doi:10.1093/oxfordjournals.jhered.a106562

CrossRef Full Text | Google Scholar

Huang, A., Xu, S., and Cai, X. (2014). Whole-genome quantitative trait locus mapping reveals major role of epistasis on yield of rice. PLoS One 9, e87330. doi:10.1371/journal.pone.0087330

PubMed Abstract | CrossRef Full Text | Google Scholar

Id-Lahoucine, S., Cánovas, A., Jaton, C., Miglior, F., Fonseca, P. A. S., Miller, S., et al. (2019a). Implementation of bayesian methods to identify SNP and haplotype regions with transmission ratio distortion across the whole genome: TRDscan v.1.0. J. Dairy Sci. 102, 3175–3188. doi:10.3168/jds.2018-15296

PubMed Abstract | CrossRef Full Text | Google Scholar

Id-Lahoucine, S., Casellas, J., Fonseca, P. A. S., Suárez-Vega, A., Schenkel, F. S., and Cánovas, A. (2022). Deviations from mendelian inheritance on bovine X-chromosome revealing recombination, sex-of-offspring effects and fertility-related candidate genes. Genes 13, 2322. doi:10.3390/genes13122322

PubMed Abstract | CrossRef Full Text | Google Scholar

Id-Lahoucine, S., and Casellas, J. (2017). Impact of incomplete pedigree data and independent culling level pre-selection on the genetic evaluation of livestock: A simulation study on lamb growth. Livest. Sci. 198, 76–81. doi:10.1016/j.livsci.2017.02.011

CrossRef Full Text | Google Scholar

Id-Lahoucine, S. (2020). Mendelian inheritance in the genomics and big data era: Transmission ratio distortion phenomenon in cattle breeds. [dissertation]. Guelph, Canada: University of Guelph. Available at: https://hdl.handle.net/10214/18142.

Google Scholar

Id-Lahoucine, S., Molina, A., Cánovas, A., and Casellas, J. (2019b). Screening for epistatic selection signatures: A simulation study. Sci. Rep. 9, 1026. doi:10.1038/s41598-019-38689-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeffreys, H. (1984). Theory of probability. Oxford: Clarendon Press.

Google Scholar

Kass, R. E., and Raftery, A. E. (1995). Bayes factors. J. Am. Stat. Assoc. 90, 773–795. doi:10.1080/01621459.1995.10476572

CrossRef Full Text | Google Scholar

Kempthorne, O. (1969). An introduction to genetic statistics. Ames Iowa: The Iowa University Press.

Google Scholar

Kiser, J. N., Keuter, E. M., Seabury, C. M., Neupane, M., Moraes, J. G. N., Dalton, J., et al. (2019). Validation of 46 loci associated with female fertility traits in cattle. BMC genomics 20, 576. doi:10.1186/s12864-019-5935-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Leppälä, J., Bokma, F., and Savolainen, O. (2013). Investigating incipient speciation in Arabidopsis lyrata from patterns of transmission ratio distortion. Genetics 194, 697–708. doi:10.1534/genetics.113.152561

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, A., Wang, Y., Sahana, G., Zhang, Q., Liu, L., Lund, M. S., et al. (2017). Genome-wide association studies for female fertility traits in Chinese and nordic holsteins. Sci. Rep. 7, 8487. doi:10.1038/s41598-017-09170-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, T. F. C. (2014). Epistasis and quantitative traits: Using model organisms to study gene-gene interactions. Nat. Rev. Genet. 15, 22–33. doi:10.1038/nrg3627

PubMed Abstract | CrossRef Full Text | Google Scholar

Maltecca, C., Weigel, K. A., Khatib, H., Cowan, M., and Bagnato, A. (2008). Whole-genome scan for quantitative trait loci associated with birth weight, gestation length and passive immune transfer in a Holstein x Jersey crossbred population. Anim. Genet. 40, 27–34. doi:10.1111/j.1365-2052.2008.01793.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Montagutelli, X., Turner, R., and Nadeau, J. H. (1996). Epistatic control of non-Mendelian inheritance in mouse interspecific crosses. Genetics 143, 1739–1752. doi:10.1093/genetics/143.4.1739

PubMed Abstract | CrossRef Full Text | Google Scholar

Muller, H. J. (1942). Isolating mechanisms, evolution and temperature. Biol. Symp. 6, 71–125.

Google Scholar

Müller, M. P., Rothammer, S., Seichter, D., Russ, I., Hinrichs, D., Tetens, J., et al. (2017). Genome-wide mapping of 10 calving and fertility traits in Holstein dairy cattle with special regard to chromosome 18. J. Dairy Sci. 100, 1987–2006. doi:10.3168/jds.2016-11506

PubMed Abstract | CrossRef Full Text | Google Scholar

Neel, J. V., and Schull, W. J. (1954). Human heredity. Chicago: University of Chicago Press.

Google Scholar

Niedzicka, M., Dudek, K., Fijarczyk, A., Zielinski, P., and Babik, W. (2017). Linkage map of Lissotriton newts provides insight into the genetic basis of reproductive isolation. G3 Genes, Genomes, Genet. 7, 2115–2124. doi:10.1534/g3.117.041178

CrossRef Full Text | Google Scholar

Ohta, T. (1982). Linkage disequilibrium due to random genetic drift in finite subdivided populations. Proc. Natl. Acad. Sci. USA. 79, 1940–1944. doi:10.1073/pnas.79.6.1940

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, P. C. (2008). Epistasis – The essential role of gene interactions in the structure and evolution of genetic systems. Nat. Rev. Genet. 9, 855–867. doi:10.1038/nrg2452

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, P. C. (1998). The language of gene interaction. Genetics 149, 1167–1171. doi:10.1093/genetics/149.3.1167

PubMed Abstract | CrossRef Full Text | Google Scholar

Sargolzaei, M., Chesnais, J. P., and Schenkel, F. S. (2014). A new approach for efficient genotype imputation using information from relatives. BMC Genomics 15, 478. doi:10.1186/1471-2164-15-478

PubMed Abstract | CrossRef Full Text | Google Scholar

Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and van der Linder, A. (2002). Bayesian measures of model complexity and fit. J. Roy. Stat. Soc. B 64, 583–639. doi:10.1111/1467-9868.00353

CrossRef Full Text | Google Scholar

Strange, T., Ask, B., and Nielsen, B. (2013). Genetic parameters of the piglet mortality traits stillborn, weak at birth, Starvation, Crushing, and Miscellaneous in crossbred pigs. J. Anim. Sci. 91, 1562–1569. doi:10.2527/jas.2012-5584

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Shang, J., and Liu, J. X. (2016). “An improved ant colony optimization algorithm for the detection of SNP-Interactions,” in Lecture notes in computer science. Editors D. S. Huang, K. Han, and A. Hussain (Switzerland: Springer), 21–32.

CrossRef Full Text | Google Scholar

Sun, Y., Shang, J., Liu, J. X., Li, S., and Zheng, C. H. (2017). epiACO - a method for identifying epistasis based on ant Colony optimization algorithm. BioData Min. 10, 23. doi:10.1186/s13040-017-0143-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomasen, J. R., Guldbrandtsen, B., Sorensen, P., Thomsen, B., and Lund, M. S. (2008). Quantitative trait loci affecting calving traits in Danish Holstein cattle. J. Dairy Sci. 91, 2098–2105. doi:10.3168/jds.2007-0602

PubMed Abstract | CrossRef Full Text | Google Scholar

Turelli, M., and Orr, H. A. (2000). Dominance, epistasis and the genetics of postzygotic isolation. Genetics 154, 1663–1679. doi:10.1093/genetics/154.4.1663

PubMed Abstract | CrossRef Full Text | Google Scholar

van Boven, M., and Weissing, F. J. (2001). Competition at the mouse t complex: Rare alleles are inherently favored. Theor. Popul. Biol. 60, 343–358. doi:10.1006/tpbi.2001.1551

PubMed Abstract | CrossRef Full Text | Google Scholar

VanRaden, P. M., Olson, K. M., Null, D. J., and Hutchison, J. L. (2011). Harmful recessive effects on fertility detected by absence of homozygous haplotypes. J. Dairy Sci. 94, 6153–6161. doi:10.3168/jds.2011-4624

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Kargo, M., Liu, A., Thomasen, J. R., Pan, Y., and Su, G. (2019). Genotype-by-environment interaction of fertility traits in Danish Holstein cattle using a single-step genomic reaction norm model. Heredity 123, 202–214.

PubMed Abstract | CrossRef Full Text | Google Scholar