TYPE Original Research PUBLISHED 06 April 2023 DOI 10.3389/fgene.2023.1132796 Parent-offspring genotyped trios OPEN ACCESS unravelling genomic regions with EDITED BY Vincenzo Landi, University of Bari Aldo Moro, Italy gametic and genotypic epistatic REVIEWED BY Antonio Molina, transmission bias on the cattle Universidad de Cordoba, Spain Mayra Gomez, Italian National Association of Buffalo genome Breeders, ANASB, Italy Qishan Wang, Samir Id-Lahoucine 1, Joaquim Casellas 2, Filippo Miglior 1, Zhejiang University, China *CORRESPONDENCE Flavio S. Schenkel 1 and Angela Cánovas 1* Angela Cánovas, 1 Centre for Genetic Improvement of Livestock, Department of Animal Biosciences, University of Guelph, acanovas@uoguelph.ca Guelph, ON, Canada, 2Departament de Ciència Animal i dels Aliments, Universitat Autònoma de SPECIALTY SECTION Barcelona, Barcelona, Spain This article was submitted to Livestock Genomics, a section of the journal Frontiers in Genetics Several biological mechanisms affecting the sperm and ova fertility and viability at RECEIVED 27 December 2022 ACCEPTED 20 March 2023 developmental stages of the reproductive cycle resulted in observable PUBLISHED 06 April 2023 transmission ratio distortion (i.e., deviation from Mendelian expectations). CITATION Gene-by-gene interactions (or epistasis) could also potentially cause specific Id-Lahoucine S, Casellas J, Miglior F, transmission ratio distortion patterns at different loci as unfavorable allelic Schenkel FS and Cánovas A (2023), combinations are under-represented, exhibiting deviation from Mendelian Parent-offspring genotyped trios unravelling genomic regions with proportions. Here, we aimed to detect pairs of loci with epistatic transmission gametic and genotypic epistatic ratio distortion using 283,817 parent-offspring genotyped trios (sire-dam- transmission bias on the cattle genome. offspring) of Holstein cattle. Allelic and genotypic parameterization for epistatic Front. Genet. 14:1132796. doi: 10.3389/fgene.2023.1132796 transmission ratio distortion were developed and implemented to scan the whole genome. Different epistatic transmission ratio distortion patterns were observed. COPYRIGHT © 2023 Id-Lahoucine, Casellas, Miglior, Using genotypic models, 7, 19 and 6 pairs of genomic regions were found with Schenkel and Cánovas. This is an open- decisive evidence with additive-by-additive, additive-by-dominance/dominance- access article distributed under the terms by-additive and dominance-by-dominance effects, respectively. Using the allelic of the Creative Commons Attribution License (CC BY). The use, distribution or transmission ratio distortion model, more insight was gained in understanding the reproduction in other forums is penetrance of single-locus distortions, revealing 17 pairs of SNPs. Scanning for the permitted, provided the original author(s) depletion of individuals carrying pairs of homozygous genotypes for unlinked loci, and the copyright owner(s) are credited and that the original publication in this revealed 56 pairs of SNPs with recessive epistatic transmission ratio distortion journal is cited, in accordance with patterns. The maximum number of expected homozygous offspring, with none of accepted academic practice. No use, them observed, was 23. Finally, in this study, we identified candidate genomic distribution or reproduction is permitted which does not comply with these terms. regions harboring epistatic interactions with potential biological implications in economically important traits, such as reproduction. KEYWORDS transmission ratio distortion, epistasis, genotyped trios, Holstein, allelic and genotypic parameterizations 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 Frontiers in Genetics 01 frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 or blocks of physically linked loci), interactions between different incomplete deviation of single-locus TRD and unraveling pairs loci (i.e., epistasis) could also potentially cause specific TRD patterns of unlinked genomic regions across the whole genome that are (Corbett-Detig et al., 2013; Leppälä et al., 2013). Within a TRD not transmitted according to Mendelian inheritance rules, but context, compatible allele (or genotype) combinations at different display epistatic transmission ratio distortion patterns. 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 2 Materials and methods addition, as epistasis influences the expressivity of the involved loci, it has been hypothesized that the incomplete penetrance of lethal 2.1 Genotypes and trios mutations may involve epistatic interaction (Ballinger and Noor, 2018). The dataset used to investigate epistatic TRD consisted of There are documented cases in model organisms, plants and fish 436,651 Holstein genotypes provided by the Canadian Dairy where epistatic interactions explain the variation of the observed TRD Network database (Lactanet, Guelph, Ontario, Canada). The number (Montagutelli et al., 1996; van Boven and Weissing, 2001; Borowsky and of genotyped trios (sire-dam-offspring) was 283,817 and were sampled Cohen, 2013; Behrouzi and Wit, 2018). Epistatic TRD at different loci from all available Holstein genotypes (>1 million; October 2017) with has been related to different sources, such as the lethality of embryos offspring genotyped within 90 days of the birth, thus minimizing carrying particular combinations of alleles (Montagutelli et al., 1996) selection artifacts (Id-Lahoucine and Casellas, 2017; Id-Lahoucine and lower fertility or diseases (Behrouzi and Wit, 2018), among others. et al., 2019a). The number of genotyped sire and dams was Particularly, epistatic TRD has been related to genic incompatibilities 5,224 and 117,316, respectively. Animals were genotyped with (or hybrid incompatibility) in crosses where deleterious epistatic different SNP genotyping arrays ranging from 6,909 to interactions between heterospecific alleles lead to hybrid sterility 777,962 SNPs (Supplementary Material S1) and mapped to the and/or inviability (Turelli and Orr, 2000; Fishman et al., 2008; UMD3.1 Bos taurus genome assembly. This data was previously Cutter, 2012; Corbett- Detig et al., 2013; Giesbers et al., 2019; imputed using FImpute (Sargolzaei et al., 2014) to 47,910 SNPs (Id- Arends et al., 2022). These genic incompatibilities between diverged Lahoucine, 2020). genomes, termed Bateson-Dobzhansky-Muller (BDM) incompatibilities, were first described by Bateson (1909), Dobzhansky (1937) and Muller (1942), representing a classic model 2.2 Analytical models statistical analyses for the evolution of reproductive isolation in diverging lineages. Despite the additive paradigm being widely predominant for the Epistatic transmission ratio distortion was evaluated on pairs of technical statistical and computational limitations, the importance SNPs across the whole genome (47,910 SNPs) using three different of interaction components of the genetic architecture of traits has methods. In order to avoid confounding due to linkage recently been emphasized (Carlborg and Haley, 2004; Phillips, disequilibrium (LD) within chromosomes and to reduce the 2008). It must be highlighted that the relative importance of computational time, epistatic TRD scan was restricted to inter- gene-by-gene interaction has been around for more than a chromosomal pairs. century (Gilbert-Diamond and Moore, 2011) and was recognized Genotypic parameterization of epistatic transmission ratio as an explanation of the observed deviations from Mendelian ratios distortion. Following Casellas et al. (2012, 2020) genotypic TRD (e.g., Bateson, 1909; Neel and Schull, 1954; Hollander, 1955; Phillips, model of one single-locus and Kempthorne (1969) epistatic theory, 1998). Identification of epistatic interaction in the genome will epistatic TRD parameters were included in the probability of the provide an opportunity to understand the epistatic effect on offspring (Poff) for each combination of genotypes of a given pair of agronomically important complex traits allowing insights into the two loci (A and B) as: genetic background and complexity underlying reproduction. The TRD approach does not require phenotypic records, but only Pof f (AABB)  1 + 1 αA + 1 αB − 1 δA − 1 δB + 1 ααe + 1 αδe + 1 δαe genomic information in trios (sire-dam-offspring), which is − 1 δδe available nowadays for most of the livestock species what opens alternative genomic strategies with potential outcomes to improve Pof f (AaBB)  2 + 0 αA + 1 αB + 1 δA − 1 δB + 0 ααe + 0 αδe − 1 δαe reproductive success. + 1 δδe The conventional method of identifying epistatic TRD relies on performing a screen for pairs of loci showing deviations from Pof f (aaBB)  1 − 1 αA + 1 αB − 1 δA − 1 δB − 1 ααe − 1 αδe + 1 δαe the expected assuming random segregation and the majority of − 1 δδe methods are restricted to an F2 design. The allelic and genotypic Pof f (AABb)  2 + 1 αA + 0 αB − 1 δA + 1 δB + 0 ααe − 1 αδe + 0 δαe TRD models developed by Casellas et al. (2012, 2014, 2017) for single-locus, which are based on tracing allele inheritance from + 1 δδe parents to offspring, are flexible and accommodate all types of Pof f (AaBb)  4 + 0 αA + 0 αB + 1 δA + 1 δB + 0 ααe + 0 αδe + 0 δαe matings and pedigree structure. Here we adapted these models to − 1 δδe account for epistatic interactions and to be applicable for livestock populations. The main objectives of this research, Pof f (aaBb)  2 − 1 αA + 0 αB − 1 δA + 1 δB + 0 ααe + 1 αδe + 0 δαe using parent-offspring genotyped trios (sire-dam-offspring) of + 1 δδe Holstein cattle, were to investigate the genetic basis for Frontiers in Genetics 02 frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 Pof f (AAbb)  1 + 1 αA − 1 αB − 1 δA − 1 δB − 1 ααe + 1 αδe − 1 δαe Kass and Raftery, 1995). More accurate estimates were obtained − 1 δδe using 550,000 iterations for the selected relevant epistatic regions from the preliminary analyses. In addition, in order to focus mainly Pof f (Aabb)  2 + 0 αA − 1 αB + 1 δA − 1 δB + 0 ααe + 0 αδe + 1 δαe on epistasis phenomenon, the ratio of maximum (BFαα, BFαδ, BFδα, + 1 δδe 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 Pof f (aabb)  1 − 1 αA − 1 αB − 1 δA − 1 δB + 1 ααe − 1 αδe − 1 δαe direct effects. The BF was also used for the multiple test correction − 1 δδe where the top 0.1% significant regions were selected as the most relevant epistatic TRD regions. where; αA and δA are additive- and dominance-TRD parameters for Allelic parameterization of epistatic transmission ratio locus A, respectively, αB and δB are additive- and dominance-TRD distortion. An allelic epistatic TRD could be targeted by tracing parameters for locus B, respectively, ααe, αδe, δαe and δδe are back the inheritance of a combination of two specific alleles of the additive-by-additive, additive-by-dominance, dominance-by- two loci (i.e., an artificial haplotype) from parents to offspring. Thus, additive and dominance-by-dominance epistatic TRD parameters, the probability of an offspring (Poff) for each combination of respectively. The parametric space considered for all the parameters genotypes of a given par of two loci (A and B) was ranges from −1 to 1. Due to the complexity of this model and the parameterized to include both direct and epistatic allelic TRD number of parameters involved, negative probabilities must be effects as: rescaled to 0 for simplification. An additional restriction is needed to guarantee Poff (AABB) + Poff (AaBB) + Poff (aaBB) + Pof f (AABB)  0.5 + βA  × 0.5 + βB  × 0.5 + βA  × 0.5 + βB  Poff (AABb) + Poff (AaBb) + Poff (aaBb) + Poff (AAbb) + Poff (Aabb) + × 1 + βAB/ab  × 1 + βAB/ab  Poff (aabb) = 1. On the other hand, it must be emphasized that all the equations of (offspring) probabilities described here correspond to Pof f (AAbb)  0.5 + βA  × 0.5 − βB  × 0.5 + βA  × 0.5 − βB  the unique case of heterozygous-by-heterozygous mating for both × 1 + βAb/aB  × 1 + βAb/aB  loci (i.e., AaBb × AaBb mating), and the probabilities needs to be Pof f (aaBB)  0.5 − βA  × 0.5 + βB  × 0.5 − βA  × 0.5 + βB  adapted for each specific mating accordingly. For instance, the Poff × 1 − βAb/aB  × 1 − βAb/aB  (AABB) becomes 4 + 1 αA + 1 αB - 1 δA - 1 δB + 1 ααe + 1 αδe + 1 δαe - Pof f (aabb)  0.5 − βA  × 0.5 − βB  × 0.5 − βA  × 0.5 − βB  1 δδe whereas Poff (AABb) is 0 for AaBB × AaBB mating. The total number of possible informative matings is 65 and which could be × 1 − βAB/ab  × 1 − βAB/ab  combined in 27 types of matings summarized in the Supplementary Poff (AABb, AAbB)  0.5 + βA  × 0.5 + βB  × 0.5 + βA  × 0.5 − βB  Material S2. × 1 + βAB/ab  × 1 + βAb/aB  Under a Bayesian implementation, the conditional posterior Poff (aaBb, aabB)  0.5 − βA  × 0.5 + βB  × 0.5 − βA  × 0.5 − βB  probabilities of the TRD parameters are defined as: × 1 − βAB/ab  × 1 − βAb/aB    pαA , δA , αB , δB , ααe , αδe , δαe , δδe y) ∝ p(y αA , δA , αB , δB , ααe , αδe , δαe , δδe  Pof f (AaBB, aABB)  0.5 + βA  × 0.5 + βB  × 0.5 − βA  × 0.5 + βB  p(αA ) p(δA ) p(αB ) p(δB ) p(ααe ) p(αδe ) p(δαe ) p(δδe ) × 1 + βAB/ab  × 1 − βAb/aB  Pof f (Aabb, aAbb)  0.5 + βA  × 0.5 − βB  × 0.5 − βA  × 0.5 − βB  where: y is the column vector of genotypes of the offspring × 1 − βAB/ab  × 1 + βAb/aB  generation. Poff (AaBb, aAbB)  0.5 + βA  × 0.5 + βB  × 0.5 − βA  × 0.5 − βB  In order to reduce the computational time, only pairs of SNPs × 1 + βAB/ab  × 1 − βAB/ab  with ≥1,000 informative offspring and ≥20 heterozygous sires or/ Poff (AabB, aABb)  0.5 + βA  × 0.5 + βB  × 0.5 − βA  × 0.5 − βB  and ≥50 heterozygous dams were used. Secondly, the observed and × 1 + βAb/aB  × 1 − βAb/aB  expected number of offspring for the 27 combined matings were determined for each pair of loci. The analyses included only pairs of where; βA and βB are direct TRD parameters for locus A and B, SNPs with ≥1,000 under- or over-represented offspring for a specific respectively, βAB/Ab, βAB/aB, βAB/ab, βAb/aB, βAb/ab and βaB/ab are the two-locus genotype offspring and mating and this difference 6 heterozygous pairwise combinations for epistatic TRD parameters. (between the observed and expected) being relatively 10% higher/ Note that AB, Ab, aB and ab are the 4 possible artificial haplotype lower than the expected number of offspring. Notice that these alleles of the two implicated SNPs. For all parameters, flat priors criteria and cut-off values were chosen based in the results observed were assumed within a parametric space ranging from −0.5 to 0.5 for in single-locus TRD analyses (Id-Lahoucine, 2020) and all these direct TRD effects and from −1.0 to 1.0 for epistatic TRD effects. The steps guaranteed a relevant deviation from Mendelian proportions parametric space of direct effect is based on the principles of in the corresponding regions allowing to target regions with Mendelian inheritance (i.e., the probability of transmission of one potential epistatic TRD signals and reducing the computational specific allele ranges from 0 (e.g., βA = −0.5) to 1 (e.g., βA = 0.5), and time considerably. The model was implemented in a Bayesian where 0.5 (e.g., βA = 0.0) corresponds to null TRD). In contrast, the approach with the metropolis-Hastings (Hastings, 1970) sampling range from −1 to 1 for epistatic parameters could be viewed as the technique. A preliminary scan was implemented with a unique positive/negative preferential transmission of one combination of Monte Carlo Markov chain of 11,000 iterations (the first two alleles (artificial haplotype) or its opposite heterozygous 1,000 were discarded for TRD estimation). The statistical combinations, respectively, whereas 0 indicates null epistatic significance of TRD parameters was tested by a Bayes factor (BF, TRD. Also, these probabilities corresponded to the unique case of Frontiers in Genetics 03 frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 heterozygous-by-heterozygous mating (i.e., AaBb × AaBb), and the 3 Results probabilities must be adapted for each specific mating accordingly. Under a Bayesian implementation, the conditional posterior 3.1 Genotypic epistatic transmission ratio probabilities of the TRD parameters are defined as: distortion   pβA , βB , βAB/Ab , βAB/aB , βAB/ab , βAb/aB , βAb/ab , βaB/aB y) ∝ p(y  βA , βB , βAB/Ab , βAB/aB , βAB/ab , βAb/aB , βAb/ab , βaB/aB ) Pairs of regions with epistatic TRD estimates were found widely pβA  pβB  pβAB/Ab  pβAB/aB  pβAB/ab  pβAb/aB  pβAb/ab  pβaB/aB  across the Holstein genome. Using a threshold of ≥1,000 for BF, a large number of pairs of SNPs were observed with decisive evidence where; y is the column vector of genotypes of the offspring for epistatic TRD. In order to discard false TRD and to identify the generation. most relevant regions, the following steps were considered. Firstly, in This method is similar to the heterozygous pairwise combination order to focus on epistasis phenomenon, regions with more procedure of haplotype analysis (single locus) described by Id- statistical relevance for direct effects than epistatic effects were Lahoucine et al. (2019a), which is highly computationally discarded. Thus, only region with a minimal ratio of demanding. For this reason, a simplified algorithm was 1,000 between the maximum BF of epistatic and direct effects implemented for a preliminary scan. The biallelic-haplotype were considered (i.e., maximum (BFαα, BFαδ, BFδα, BFδδ)/ procedure for haplotype analysis described by Id-Lahoucine et al. maximum (BFαA, BFδA, BFαB, BFδB) > 1,000). Secondly, estimates (2019a) was adapted to estimate TRD for the artificial haplotypes in with large credible intervals were discarded as artifacts of the order to target signals of epistasis. Thus, the transmission probability convergence of the model following Id-Lahoucine (2020). For (P) of an artificial haplotype from heterozygote parents to offspring this purpose, pairs of SNPs with a coefficient of variation >20% was parameterized including one overall epistatic TRD effect (αj) for for significant epistatic effects were excluded. Notice that from each specific j artificial haplotype (AH) allele (P (AHj) = 1—P (AH-j) = single-locus TRD, when clear TRD exist, the standard variation 0.5 + αj). of TRD estimates is mostly null (<0.01), thus, the use of coefficient of As an exhaustive search for epistatic allelic TRD is variation of TRD estimates is a straightforward rule to discard computationally unfeasible, epistatic overall TRD were estimated regions with instable convergence. Thirdly, given that most only for the combination of pairs of SNPs which include one SNP regions have several direct and epistatic effects with different with significant direct TRD effect already identified from analysis of statistical significance, the results were separated into their single-locus TRD (2,962 SNPs; Id-Lahoucine, 2020). After the corresponding effects (i.e., additive-by-additive, additive-by- construction of artificial bi-allelic haplotypes, the analyses were dominance/dominance-by-additive or dominance-by-dominance) performed within a Bayesian framework using a TRDscan according to the most relevant effect in terms of BF. After v.1.0 software (Id-Lahoucine et al., 2019a) with a unique Monte filtering the results following the previous steps, the number of Carlo Markov chain of 110,000 iterations where the first the obtained pairs were 59,831, 87,699 and 20,549 for additive-by- 10,000 iterations were discarded as burn-in. The statistical additive, additive-by-dominance/dominance-by-additive and significance of TRD was evaluated using a Bayes factor (Kass and dominance-by-dominance effects, respectively. Raftery, 1995). Following this analysis, the candidate regions with After implementing a multiple test correction, which was based signals of epistatic allelic pattern obtained by the simplified method on selecting the top 0.1% of the regions within each category were re-examined with the full epistatic allelic model (with according to BF, the number of regions reduced to 169, being 60, 550,000 iterations), thus, allowing to overcome the limitation of 88 and 21 with additive-by-additive, additive-by-dominance/ this simplified method which was not parameterized to take into dominance-by-additive and dominance-by-dominance effects, account both direct and epistatic effects simultaneously. Finally, respectively. This strategy was implemented with the objective to both allelic and genotypic parameterizations were compared using select the most significant “top hits” of epistatic TRD signals the deviance information criterion (DIC, Spiegelhalter et al., 2002) following François et al. (2016) and Boschiero et al. (2018). to determine the goodness-of-fit and the epistatic pattern of each These pairs of regions were chosen for a second analysis with a pair of regions. large Monte Carlo Markov chain of 500,000 iterations. Applying the Double recessive epistatic transmission ratio distortion. The previous criteria described above to the new estimated, 55 out of criteria applied for both genotypic and allelic model will discard the 88 regions with additive-by-dominance/dominance-by-additive possibility of detecting double recessive epistatic TRD. In this sense, were discarded. Moreover, among the identified pairs, we in order to target double recessive epistatic TRD, the biallelic- observed that some SNPs interacted with more than one haplotype procedure with artificial haplotypes was implemented individual SNP. For instance, the SNP BTA17: to screen for absent of double homozygous offspring for two-locus 39,696,262 interacted with 26 SNPs covering from 4,197,354 bp genotypes. For this, a minimal 15 none-observed offspring for one to 12,447,484 bp on BTA23. Within this context, only the pair of single combination of homozygous genotypes were chosen, while SNPs with the highest BF were maintained as the best candidate ensuring no deviation from Mendelian proportions for other regions when one single locus interacts with several physically linked genotypes. SNPs. Thus, the number of regions reduced to 7, 19 and 6 with Multiple test correction: For multiple test correction, the BF additive-by-additive, additive-by-dominance/dominance-by- was used in order to choose the best candidate regions with additive and dominance-by-dominance effects, respectively. These epistatic effects. Thus, the multiple test correction was performed results are summarized in Table 1 for additive-by-additive and by selecting the top 0.1% of the regions within each category dominance-by-dominance effects and in Table 2 for additive-by- according to BF. dominance/dominance-by-additive. Additional details (e.g., Frontiers in Genetics 04 frontiersin.org Frontiers in Genetics Id-Lahoucine et al. 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. Coordinates locus αAa αB (Log10 δA δB ααe αδe δαe δδe Log10 DICd genotypic DIC allelic A and B (Log10(BFb)) (BF)) (Log10 (Log10 (Log10 (Log10(BF)) (Log10(BF)) (Log10(BF)) (LRc) model model (BF)) (BF)) (BF)) 1:41490266 × 21: 0.15 (693.7) −0.04 (22.8) 0.43 0.44 1 (6794.4) −0.22 (850) −0.03 (9.9) −0.94 (6459.6) 12033 269451.4 283614.8 26405185 (1524.4) (1756.7) 1:95213011 × 23: −0.06 (41.4) −0.29 0.44 0.33 (853.2) −1 (5112.1) −0.07 (46.3) 0.42 (1663.7) −0.88 (4899.6) 10342.8 243497.0 265045.6 30121507 (1358.8) (2169.2) 4:85759993 × 21: 0.38 (1191.4) −0.24 0.22 (366.2) 0.29 (639.9) 1 (4681.5) −0.43 (1498.1) 0.25 (736.4) −0.64 (2876.6) 7971.1 272790.3 278502.2 37978844 (776.9) 10:20146107 × 14: 0.3 (1014.1) −0.06 0.33 (915.8) 0.42 1 (5903.5) −0.4 (1748.2) −0.01 (−0.1) −0.82 (4587) 10775.9 266355.8 278134.6 29246158 (397.4) (1846.2) 17:39696262 × 23: 0.05 (24.5) −0.27 0.46 (3596) 0.4 (1365.5) 1 (6648.6) 0.06 (73.2) 0.37 (1619.3) −0.94 (6250.2) 12782.1 268513.8 283693.9 5524345 (1239.7) 21:60175026 × 27: 0 (−1.5) 0.18 (896.7) 0.34 (862.5) 0.3 (670.5) 1 (4880.7) −0.05 (23.6) −0.23 (889.3) −0.74 (3597.7) 8058.9 264658.6 276906.1 1343227 23:1674622 × 28: 0.57 (605.3) −0.67 −0.36 (450) −0.6 −0.99 0.51 (598.9) −0.44 (410.9) −0.18 (114.5) 5117 215691.1 262136.3 05 2233574 (759.9) (1355.6) (5388.6) 1:36209316 × 21: 0.26 (1526.4) 0.14 (756.9) 0.48 0.47 1 (10097.7) −0.34 (2110.8) −0.24 (1429.7) −1 (16219.5) 25721.6 235138.9 272410.9 26405185 (1906.4) (9303.7) 1:95101691 × 23: 0.09 (66.9) −0.06 (23.7) 0.39 0.33 (979.3) 1 (5376.7) 0.11 (111) 0.21 (829.7) −1 (13060.3) 15124.9 239767.8 268110.2 29897088 (1362.2) 4:85384769 × 21: 0.45 (3133.8) 0.33 0.4 (4817.4) 0.47 1 (9437.1) −0.54 (3956.6) −0.49 (3653.7) −1 (20223.3) 32514.9 202493.3 252222.1 37978844 (2055.7) (4527.9) 10:25732248 × 14: −0.18 (863.3) 0.44 (2600) 0.58 0.39 (9252) −1 (9299.1) 0.33 (1901.5) −0.62 (4533.1) −1 (15595.4) 27517.3 214749.7 257410.0 29246158 (6551.3) 13:1655502 × 28: −0.1 (49.7) 0.17 (159.4) 0.14 (148.6) 0.19 (240.6) 0.8 (1708.3) 0.13 (96.1) −0.11 (76.8) −0.59 (8447) 3613.9 267493.0 279825.4 21945694 17:39696262 × 23: 0.2 (950) −0.37 0.51 0.38 1 (7685.7) −0.1 (75.5) 0.52 (2916.7) −1 (10688.5) 18451.2 247479.4 276381.9 10.3389/fgene.2023.1132796 6727000 (1855.9) (8181.6) (1435.3) a α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. b Bayes Factor. c Likelihood ratio [= L (y | αA, δA, αB, δB, ααe, αδe, δαe, δδe)/L (y | αA = 0, δA = 0, αB = 0, δB = 0, ααe = 0, αδe = 0, δαe = 0, δδe = 0)]. frontiersin.org d Deviance information criterion. Frontiers in Genetics Id-Lahoucine et al. TABLE 2 Pairs of genomic regions identified with epistatic transmission ratio distortion (TRD) in Holstein cattle with mainly additive-by-dominance effects. Coordinates locus αAa αB (Log10 δA δB ααe αδe δαe δδe Log10 DICd genotypic DIC allelic A and B (Log10 (BF)) (Log10(BF)) (Log10(BF)) (Log10(BF)) (Log10(BF)) (Log10(BF)) (Log10(BF)) (LRc) model model (BFb)) 1:33593708 × 21:26405185 0.43 (357) −0.71 0.11 (751.5) −0.12 (670.2) 1 (1443.4) −0.73 (1052.1) 0.92 (1917.7) −0.52 (1220.9) 5319.8 221678.0 219188.4 (1273.1) 1:44623201 × 7:10835967 0.72 (3173.9) 0.02 (0) −0.15 (184.7) 0.64 (2276.9) 0.03 (3.5) −0.47 (4339.9) 0.02 (0.8) −0.07 (36.9) 3628.4 246495.6 264314.5 2:27724930 × 7:10835967 0.4 (1073.9) 0.03 (2.1) −0.04 (8.1) 0.67 (1566.7) 0.04 (3.6) −0.52 (1646.3) −0.01 (−1.4) 0 (−1.3) 2753.3 209389.9 191723.9 4:18207858 × 7:10835967 −0.52 0.08 (44.2) −0.07 (38.7) 0.7 (2655.6) −0.02 (0.5) 0.36 (2845.3) −0.03 (4.4) −0.1 (82.4) 3372 247299.9 261595.1 (1757.5) 4:85188654 × 21:37978844 0.23 (42.4) 0.67 (1067.4) 0.57 (1730.8) 0.12 (89.5) −1 (1473) −0.27 (122.7) −0.93 (2212.7) −0.77 (2208.7) 5816.1 176258.6 178503.0 5:2793760 × 17:72942592 −0.9 (2033.2) 0.15 (88.8) −0.2 (276.7) 0.29 (359.8) −0.21 (128.9) 0.3 (2082.9) 0.03 (1.4) −0.01 (−1.1) 2352.4 237918.0 250837.9 5:56169416 × 7:10835967 −0.5 (1508.3) 0.06 (12.1) 0 (−1.4) 0.63 (1429) −0.06 (7.6) 0.56 (1931.3) −0.02 (−0.4) −0.06 (19.8) 2519.5 220619.4 212445.1 6:4076731 × 7:10835967 0.52 (1493.3) 0.06 (23.6) −0.09 (66.8) 0.65 (2364.1) 0.03 (3.4) −0.35 (3647.7) −0.02 (0.4) −0.08 (51.6) 3048.1 238330.7 258680.7 7:10835967 × 8:101251865 0.02 (0.8) −0.46 (787.9) 0.69 (1468.1) 0 (−1.1) −0.02 (−0.2) 0 (−1.4) 0.62 (1545.9) −0.04 (6) 2833.8 199930.3 179384.2 7:10835967 × 10:29761646 0.07 (29.3) 0.69 (3299.7) 0.67 (2574.4) −0.19 (361.4) 0.04 (7.2) −0.01 (−1.5) −0.44 (3652.7) −0.04 (14.5) 3957.7 256869.6 272963.3 06 7:10835967 × 11:77484457 0.06 (22.9) 0.58 (2209) 0.64 (2631.8) −0.06 (27.9) 0.03 (3.1) −0.02 (0.1) −0.38 (4236.9) −0.11 (117) 3438.7 262239.4 281357.1 7:10835967 × 16:34911228 0 (−1.4) 0.38 (889.3) 0.68 (1568.8) −0.02 (−0.1) 0.04 (5.9) 0 (−1.4) −0.57 (1609.4) −0.05 (16.3) 2815.2 207906.1 187295.0 7:10835967 × 17:48671784 0.05 (10.8) 0.42 (1673.5) 0.63 (1491.3) −0.07 (24.1) 0.03 (0.8) −0.01 (−1) −0.53 (1737.3) 0.01 (−1.3) 2669.2 218822.9 207813.4 7:10835967 × 21:44499384 0.06 (15.4) −0.37 0.66 (1537.6) 0 (−1.3) −0.01 (−0.6) −0.02 (1.1) 0.51 (1486.6) −0.03 (5.1) 2644.1 211959.4 192902.4 (1084.4) 7:10835967 × 24:33863680 0.02 (0.3) −0.37 0.68 (1673) −0.01 (−0.3) −0.03 (2.6) −0.01 (−1) 0.51 (1559.1) −0.05 (11.8) 2755.5 210855.8 194410.1 (1189.2) 8:76620508 × 25:28413439 −0.47 (44.8) −0.97 (477.8) −0.15 (16.7) 0.43 (180.9) 0.06 (4.1) 0.74 (1663.6) −0.18 (20.3) 0.3 (99) 551.9 112941.5 105642.5 9:59595970 × 18:28139889 −0.57 (36.8) −0.78 (313.5) −0.31 (35.2) 0.48 (212.2) 0.39 (114.6) 0.84 (1223.8) −0.35 (75.2) 0.5 (246.6) 247.6 122457.6 111722.0 10:26700563 × 14:29246158 −0.49 (455.7) −0.62 (898.8) −0.03 (13.5) −0.1 (550.8) −1 (1667.1) 0.85 (1627.5) 0.88 (1964.1) −0.5 (1132.9) 5973.3 235802.2 234589.0 15:82260685 × 18:10642175 −0.55 0.18 (105) 0.04 (8.6) 0.04 (6.1) 0.72 (1109.4) 0.23 (305.8) −0.33 (2251.6) −0.25 (608.7) 2261.2 258076.8 270979.5 10.3389/fgene.2023.1132796 (1299.8) a α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. b Bayes Factor. Likelihood ratio (= L (y | αA, δA, αB, δB, ααe, αδe, δαe, δδe)/L (y | αA = 0, δA = 0, αB = 0, δB = 0, ααe = 0, αδe = 0, δαe = 0, δδe = 0)). frontiersin.org c d Deviance information criterion. Id-Lahoucine et al. 10.3389/fgene.2023.1132796 number of informative parents, frequencies, etc.) are presented in TRD effects, displayed significant epistatic TRD effect with the Supplementary Material S3. log10(BF) > 9.61, being log10(BF) = 969.06 the maximum value When comparing between the genotypic and allelic models, observed for an epistatic heterozygous pairwise combination. The different goodness-of-fit values were observed among the regions number of regions that showed clear epistatic effects with the full with additive-by-additive effects displaying reductions from model in all heterozygous pairwise combinations including the 1,186.35 to 46,445.15 DIC units favoring the genotypic model artificial haplotype were 29 and reduced to 17 when only one over the allelic model. Notice that models with smaller DIC single artificial haplotype displayed epistatic TRD (Table 3). 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 3.3 Recessive epistatic transmission ratio dominance-by-dominance effects, reductions up to 49,728.74 DIC distortion units were observed. In addition, out of the regions found with additive-by-dominance/dominance-by-additive, only 8 pairs of The total number of pairs of regions with at least 15 expected regions favored the genotypic model with differences of up homozygous offspring but none of them observed assuming random 20,349.96 DIC units. Therefore, the remaining pairs of regions assortment for both implicated SNPs, was 67. This was reduced to (11) detected with the genotypic model with additive-by- 56 after discarding pairs of regions pointing to physically related dominance/dominance-by-additive showed better goodness-of-fit regions (Table 4, Supplementary Material S5). The maximum with the allelic model in terms of DIC. Nevertheless, more number of expected homozygous offspring for both SNPs, but statistical relevance based on BF was observed for TRD none of them observed, was 23. The number of informative sires, parameters of genotypic model (up to log10(BF) = 1,964.13) in dams and offspring for the identified regions reached to 105, comparison to TRD estimated of allelic model (up to log10(BF) = 2,809 and 11,995, respectively. 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., 4 Discussion additive-by-additive, additive-by-dominance) and BF. The current accessibility of high-throughput genotyping technologies in big data era have allowed the investigation of 3.2 Allelic epistatic transmission ratio epistasis more deeply. Previous studies used genotype data and distortion phenotypes or gene expression to investigate epistatic interactions between loci (e.g., Strange et al., 2013; Hemani et al., 2014; Huang Using the simplified method of the allelic model, a minimal et al., 2014; Mackay, 2014). Nevertheless, though many efforts have increase of 0.05 in the TRD magnitude explained by the artificial been done for identifying epistatic interactions, the methods/ haplotype (i.e., gamete) compared to the single SNP and an equal or algorithmics are still being developed (e.g., Beissinger et al., 2015; greater number of under-/over-represented offspring (i.e., ≈ αj × 2 × Sun et al., 2016; 2017; Behrouzi and Wit, 2018; Id-Lahoucine et al., number of informative offspring) were considered to identify 2019b). In this paper we extended TRD models from single-locus to candidate regions with allelic epistatic TRD pattern. The total capture epistatic interactions, presenting a practical methodology number of obtained regions detected with BF ≥ 1,000 was 4,852. when trios of parent-offspring genotypes (sire-dam-offspring) are For the identified artificial haplotypes, the absolute TRD magnitude available, not being restricted to only F2 designs, which is the most ranged from |0.06| to |0.42|, and the maximum number of under- commonly used design in recent studies of TRD and epistasis in and/or over-represented offspring was 12,756 genotypes. The model organisms (e.g., Brennan et al., 2014; Niedzicka et al., 2017; additional magnitude of TRD explained by artificial haplotypes Haddad et al., 2018). 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: 4.1 Genotypic epistatic transmission ratio 6,948,746) displaying an overall TRD of −0.22 and 3,078 under- distortion represented offspring (log10(BF) = 302.62) and exhibited a deviation of −0.37 and 3,908 under-represented offspring (log10(BF) = 692.38) The Bayes factor is the standard Bayesian tool to compare two when paired with a specific allele in another SNP (BTA1: competing models (Kass and Raftery, 1995). For determining 45,157,959). Thus, when the same SNP was found interacting in statistical significance, the possible choices of a BF threshold several pairs, we considered the pair of loci with the highest BF as the could be variable, BF ≥ 10 indicating strong evidence and BF ≥ best candidate region with epistatic TRD. This reduced the number 100 more decisive evidence according to the Jeffreys’ (1984) scale. of pairs of regions to 78 (Supplementary Material S4). In terms of Within the context of TRD analysis, the BF can be viewed as a goodness-of-fit, among the 78 pairs of regions detected with the measure of the strength of statistical evidence of the detected allelic model, minimal DIC values were observed on 76 pairs when regions, given both TRD magnitude and sample size of compared to the genotypic model. Specifically, reduction from informative offspring simultaneously. The value of BF increases 19.38 up to 23,647.43 DIC units were observed for the allelic with the increase of TRD magnitude and the number of informative model in comparison to the genotypic model. The re-analyses of offspring (Id-Lahoucine et al., 2019a). Here, a higher threshold of BF these regions with the full allelic model, with epistatic and direct (≥1,000) was considered to ensure targeting epistatic regions with Frontiers in Genetics 07 frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 relevant TRD magnitudes. Even though the prevalence of signals of showed an under-representation of 14,338, 14,079, 25,459 and epistatic TRD widely extended across the Holstein genome (i.e., a 24,814 for AaBB, AABb, aaBb and Aabb, respectively, and an large number of pairs of SNPs (>150,000) were initially identified), it over-representation of 17,864, 1, 37,346, 9 and 29,853 for AABB, is expected that part of them are simply artifacts of the sampling aaBB, AaBb, AAbb and aabb, respectively, following the fluctuations that generate random/false TRD. It should be noted that parameterization for dominance-by-dominance effects (Figure 1B, the complexity and dimensionality of the model itself with Supplementary Material S3). Those examples are the simplest cases 8 parameters and their possible interactions in the same data with clear epistatic TRD patterns, but it must be taken into could generate false epistatic TRD. It was previously reported consideration that most regions have several direct and epistatic that the structure of the data with an unbalanced number of trios effects with different statistical significance as evidence of the across matings could easily generate TRD artifacts (Casellas et al., complexity of the epistatic phenomenon in the genome of cattle. 2020; Id-Lahoucine, 2020). For example, an inappropriate Furthermore, identification of SNPs interacting with several convergence of TRD estimation is produced when combinations physically linked SNPs supports the relevance of the epistasis of different TRD parameters maximize the likelihood of the data and phenomenon. This observation is potentially explained by the result in large credible intervals for TRD estimates (including often LD, where one individual or different physically linked SNPs zero values). For this reason, a minimal dispersion of TRD effects (i.e., a genomic region/locus) interact with one or different was assumed to discard potential TRD artifacts. physically linked SNPs, supporting the importance of the Given the scope of the study, which focus mainly on epistasis epistatic interaction found. For this reason, only the pair of SNPs phenomenon, the ratio of maximum (BFαα, BFαδ, BFδα, BFδδ)/ with the highest BF were maintained as the best candidate regions maximum (BFαA, BFδA, BFαB, BFδB) was used to select the when one single locus interacted with several physically linked SNPs. regions with more statistical evidence of epistatic effects rather On the other hand, an individual SNP interacting with several not direct effects. Assuming, for example, that BFαα [= p (ααe = 0)/p linked SNP were also found. This latter case might potentially (ααe = 0|y)] and BFαA [= p (αA = 0)/p (αA = 0|y)] are the maximum suggest a multiple interaction involving different unlinked loci values, this ratio becomes = p (αA = 0|y)/p (ααe = 0|y). Given p (αA = (of order 3 or more). 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 4.2 Allelic epistatic transmission ratio ensured that the epistatic effect is much more relevant to explain the distortion observed data of the implicated regions than the direct effects. After the multiple test correction (top 0.1% of the regions within The preliminary scan with the allelic model was implemented each category according to BF), the minimal and maximal log10(BF) only in pairs of loci with already identified SNP with direct TRD effect obtained after the multiple test correction was 1,223.8 and 20,223, in order to target possible epistatic interaction explaining the respectively. Notice that BF is a measure of the strength of statistical incomplete penetrance (deviation) of single-locus TRD. The evidence (for both TRD magnitude and sample size of informative criteria implemented in the simplified method were used to ensure offspring). The likelihood of the data assuming epistatic effects for more TRD was explained with epistatic interaction and that the these regions is ≥ 101223.8 times more probable that assuming null observed TRD of the individual SNP was captured by an artificial TRD effects, very strongly supporting the presence of epistatic haplotype. Explicitly, the restriction of an equal or greater number of effects. under-/over-represented offspring was used to guarantee that the It is important to mention that while the number of regions with under-/over-represented offspring generated from a specific parent additive-by-additive and dominance-by-dominance did not reduce with a particular SNP allele corresponds to a single artificial haplotype after using large sample chains, they were considerably reduced for (gamete) on those parents. Note that the linkage of specific alleles in a regions with additive-by-dominance/dominance-by-additive effects. single gamete could be a result of selection or random effects, as it is This latter observation is partially due to the confounding between the known by theory that selection or genetic drift could favor direct dominance effect and the additive-by-dominance/dominance- combinations of alleles and consequently induce LD at physically by-additive epistatic effects as observed in different estimates among unlinked loci (Bulmer, 1971; Ohta, 1982). However, by using trios of different chains in some cases. This suggests that several combinations parent-offspring genotypes with offspring genotyped at early age, the of TRD estimates maximize the likelihood of the data, as mentioned issue of selection potentially was minimized, and the detected regions before, and therefore, it is more difficulty to estimate their effects. are expected to be associated with the reproduction cycle. In this case, Moreover, regions with additive-by-additive effects showed the under the hypothesis that the partial lethality of one allele at the first clearest patterns of epistatic effects. One of the regions with an locus is linked to a specific allele on the second paired locus, this estimated ααe = 0.99 (standard deviation = 0.00015, log10(BFαα) = methodology could allow more insights to be gained in understating 4,880.73) displayed an over-representation of 8,054 offspring (40.04%; the penetrance of lethal alleles when knowing the genotype of the 20,117 observed vs. 12,063 expected) for AABB and 14,800 (30.29%; implicated loci. For instance, the variation of penetrance explained by 48,867 observed vs. 34,067 expected) for aaBB, whereas an under- the interaction could be due to some functional interaction of some representation of 5,567 (49.80%; 5,633 observed vs. 11,179 expected) unfavorable allelic combinations when inherited together in the for aaBB and 4,227 (66.85%; 2,100 observed vs. 6,323 expected) for genome of the progeny and probably affecting the fertility of AAbb (Figure 1A, Supplementary Material S3). On the other hand, parents and/or viability of embryos/offspring. regarding the dominance-by-dominance effect, a pair of regions with On the other hand, it is important to mention that, as the δδe = −0.99 (standard deviation = 0.00004, log10(BFδδ) = 16,219.5) analyses are restricted to SNPs with single-locus TRD in which the Frontiers in Genetics 08 frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 majority had low frequencies (Id-Lahoucine, 2020), the multiple test (minimal observed was 1012.71). More specifically, the epistatic combination with all other SNPs could easily generate allelic parameters seems to be very importance as the likelihood ratio epistasis in some of them by chance. Indeed, SNPs interacting between the full model and an alternative model without epistatic with several SNPs were observed. Specifically, only 11 pairs of parameters exhibited values up to 101,261.89 (minimal 2,630.27), SNPs had one single interaction whereas the maximum number supporting the importance of the epistatic TRD phenomenon in of SNPs observed interacting with one individual SNP was 866. Even the data. though it is plausible to assume that there are multiple interactions Finally, it is important to highlight that not always a between several loci, these multiple interactions must be taken with disadvantage was observed when a specific artificial haplotype caution given the low frequencies of the analyzed SNPs (from the was transmitted, but also artificial haplotype alleles with 2,962 used SNPs, 1,424 and 950 had frequency <0.05 and 0.01, positive effects were found. In this sense, one single artificial respectively). However, given the analysis were restricted to haplotype with a negative TRD was observed, associated with candidate regions (with low frequencies in most cases), the pair an under-represented offspring and potentially lethality of the of loci with the highest BF was considered as the best candidate offspring of carrier parent, whereas the alternative artificial region with epistatic TRD, when the same SNP was found haplotypes had positive effects. In contrast, when one single interacting with several SNPs. artificial haplotype had a positive effect (whereas negative or The results of the analysis using the full model including null for the remaining artificial haplotypes), it suggests the simultaneously both direct and epistatic effects supported the advantage of one haplotype to be transmitted over other previous findings. First, this emphasizes the effectiveness of the haplotypes. This could be due to an increase in fitness when an simplified approach to largely reduce the computational time. allele is paired with a specific allele at another locus. In addition, it Second, in comparison to the simplified method, the epistatic could be hypothesized that these regions with positive selection are effects observed in one artificial haplotype were also observed in more related to biological processes acting before fertilization, all (3) or some of the heterozygous pairwise combinations including where gametes with specific allele combinations have a preferential the artificial haplotype identified. An example of this was an artificial transmission to the next-generation. haplotype (aB) with a TRD magnitude of 0.37 (log10(BF) = 769.12) that, using biallelic-haplotype method, showed a TRD of −0.78 (log10(BF) = 351.96), −0.67 (log10(BF) = 216.05) and 0.74 4.3 Recessive epistatic transmission ratio (log10(BF) = 192.72) for βAB/aB, βAb/aB and βaB/ab, respectively, distortion whereas keeping βA = −0.31 (log10(BF) = 471.30) and βB = 0.29 (log10(BF) = 328.74) for direct effect. The number of regions which Absence (or depletion) of individuals with two-locus genotypes showed clear epistatic effects in all heterozygous pairwise in the homozygous statue were found across the Holstein genome, combinations including the artificial haplotype, was 29. suggesting multiple variants with recessive epistatic TRD pattern. Nevertheless, in some cases, epistatic TRD was also observed in The number of candidate pairs of loci were 56, with 23 of them heterozygous pairwise combinations with other artificial haplotype having the maximum number of expected homozygous offspring, alleles. This result adds difficulties in understanding epistatic TRD, but none of them being observed. It is important to mention that but could be explained by the different linkage with the causal Holstein haplotype 3 was identified initially by VanRaden et al. mutations among families, higher order of epistasis interactions, and (2011) with only 7 non-observed homozygous offspring from does not rule out the possibility of involvement of specific heterozygous sires in combination with heterozygous maternal interactions between the artificial haplotypes. A total of grandsires. These findings provide evidence of regions carrying 17 regions were identified when only pairs of regions displaying deleterious mutations that are expressed when presented together epistatic TRD for a single artificial haplotype were considered in a homozygous state, potentially producing lethality. (Table 3). Moreover, we observed that the magnitude obtained from single-locus analyses were similar (mostly identical) to those 4.4 Quantitative trait loci and epistatic obtained from the epistatic allelic model (direct effects) when transmission ratio distortion both direct and epistatic TRD effect were simultaneously estimated. This interesting behavior support the robustness and To support TRD findings, several quantitative trait loci (QTL) the accuracy of allelic model, where basically the model it tried to from the Cattle QTL database (CattleQTLdb, www.animalgenome. track the inheritance of alleles from parents to offspring org/cgi-bin/QTLdb/BT/index) were found to overlap with epistatic independently of its source. It is also important to mention that TRD regions. Among them, the pair of SNPs BTA1:36209316 and a similar behavior was already reported in single-locus TRD BTA21:26405185 were reported with QTL for reproduction traits analyses, thus highlighting the robustness of allelic model when such as gestation length (Maltecca et al., 2008) and first service to comparing both allelic and genotypic model (Id-Lahoucine, 2020; conception traits (Kiser et al., 2019). Cole et al. (2011) reported QTL Id-Lahoucine et al., 2022). Thus, this latter suggests that epistatic for both Calving ease and Stillbirth for BTA23:1674622 and BTA28: TRD could generate signals of single-locus TRD (as expected) and 2233574. Other pairs of SNPs (BTA 2:27724930 and BTA 7: consequently could make difficult to differentiate between the 10835967) were also found with QTL for calving ease and origins of single-locus TRD in some cases. However, it is stillbirth by Cole et al. (2011); Müller et al. (2017). In addition, important to emphasize that the likelihood ratio between the QTL for the interval from first to last insemination and non-return full and null TRD models exhibited values up to 101448.11 rate (Müller et al., 2017), interval to first estrus after calving (Liu Frontiers in Genetics 09 frontiersin.org Frontiers in Genetics Id-Lahoucine et al. TABLE 3 Pairs of genomic regions identified with epistatic transmission ratio distortion (TRD) in Holstein cattle with allelic patterns. Coordinates Coordinates Artificial overall-TRD βAc βB βAB/Ab βAB/aB βAB/ab βAb/aB βAb/ab βaB/ab DICd DIC locus A locus B haplotypea for AH (Log10(BF)) (Log10(BF)) (Log10(BF)) (Log10 (Log10(BF)) (Log10(BF)) (Log10 (Log10(BF)) allelic genotypic (Log10(BFb)) (BF)) (BF)) Model Model 1:26453228 10:29751813 aB −0.13 (54.4) 0.03 (3.7) 0 (−2) 0.01 (−1.7) 0.08 −0.31 (42.7) 0.33 (41.6) 0 (−0.9) −0.4 (0.3) 119791.6 135155.3 (1.6) 1:39339779 7:94639822 Ab −0.18 (105.4) 0 (−2) 0.05 (14.7) 0.47 (58.3) 0 (−1.7) −0.27 (35.4) −0.17 (11.6) 0 (0)4 −0.07 (1.8) 135191.8 156489.9 1:45157959 23:6948746 Ab −0.37 (692.4) 0.14 (98.3) 0.22 (302.6) 0.73 (693.5) 0.28 0 (0) 0 (0) 0 (0) 0 (0) 5093.482 5112.857 (90.2) 2:16177570 12:38409372 aB −0.27 (104.7) 0 (−2) −0.12 (38.6) −0.09 (0.3) −0.5 0.3 (13.2) 0.34 (15.8) 0 (−1.7) −0.58 (57.6) 118348.2 134441.9 (0.4) 2:116064542 19:15841374 ab −0.17 (107.7) 0.01 (4.3) 0.01 (3.7) 0 (−1.7) −0.01 0.1 (5.5) −0.14 (11.3) 0.92 0 (-0.5) 71487.01 73503.76 (−1.2) (236.3) 3:115888900 15:41707747 aB −0.09 (109.5) 0.02 (19.5) −0.01 (1.9) 0 (−1.7) 0.14 −0.14 (40.5) 0.09 (15.2) −0.01 −0.05 (-0.7) 151586.2 161556.2 (16.5) (-1) 4:32529693 13:37865729 AB −0.26 (255.6) −0.1 (62.8) 0 (−2) 0.57 (0.7) −0.57 −0.36 (50.2) 0.38 (63.9) −0.02 0 (−1.7) 134058.1 155799.7 (117.2) (−1) 5:52771740 7:23568940 ab −0.11 (81) 0 (−2) 0.04 (17.8) −0.05 (−0.3) 0 (−1.7) 0.17 (15.6) 0.25 (42.9) 0.01 (0) 0.08 (5) 134919.9 154360.6 5:73713772 12:5817779 ab −0.4 (313.7) 0.23 (149.1) 0 (−2) −0.01 (−1.7) 0.25 0.77 (216) −0.22 (8.5) 0.05 0.51 (0.4) 129106.5 149538.6 10 (4.1) (-0.3) 5:106734118 8:9600065 AB −0.17 (78.9) −0.07 (20.7) 0 (−2) 0.08 (-0.4) −0.31 −0.29 (24.4) 0.27 (21.7) 0.04 0.01 (−1.7) 110051.4 123003.1 (27.6) (−0.7) 7:76349643 11:77560928 aB −0.12 (135.7) 0.06 (47.6) 0 (−1.9) 0 (−1.7) 0.1 −0.21 (52.2) 0.23 (59.6) 0.18 0 (-0.2) 137400.3 154216.9 (10.9) (6.4) 8:33964937 21:19601568 Ab −0.12 (54) 0 (−2) 0.02 (2.7) 0.12 (3.7) 0 (−1.7) −0.21 (28.2) −0.2 (22) −0.2 −0.09 (6.8) 136738.1 156144.7 (0.1) 8:49283003 10:95891436 Ab −0.1 (39.5) 0 (−2) 0.05 (12.1) 0.12 (8.4) 0 (−1.7) −0.22 (14.8) −0.31 (28.7) −0.6 −0.01 (−0.9) 104062.9 114919.1 (0.6) 9:44932317 27:30050354 ab 0.27 (99) 0.04 (7.6) −0.01 (−0.4) 0 (−1.7) 0.06 −0.57 (107.4) 0.62 (137.9) 0.09 −0.6 (0.6) 73645.82 77756.45 (1.3) (-0.4) 9:86330663 24:9732103 Ab −0.1 (122.8) −0.01 (4.3) 0.04 (37.8) 0.3 (88.9) −0.01 −0.13 (25.8) −0.04 (1.5) 0.03 −0.04 (1.8) 148273.6 166832.4 10.3389/fgene.2023.1132796 (−1.7) (-0.6) 13:18832047 29:48033872 aB −0.24 (54.3) 0 (−2) −0.14 (30.9) 0.02 (−0.9) 0 (0) −0.04 (−0.8) 0.08 (-0.4) 0 (−1.7) −0.53 (39.2) 129725.9 149664.9 16:24948459 21:42668952 ab −0.18 (46.6) 0.07 (10.7) 0 (−2) −0.01 (−1.7) 0.07 0.43 (40.9) −0.26 (12.3) 0.13 −0.01 (0.1) 111212.7 125056.7 (−0.6) (1.5) frontiersin.org a Biallelic-haplotype method. b Bayes Factor. c β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, model. d Deviance information criterion. Frontiers in Genetics Id-Lahoucine et al. TABLE 4 Pairs of genomic regions identified with epistatic transmission ratio distortion (TRD) in Holstein cattle with recessive patterns. Coordinates Coordinates Artificial haplotype number of number of Frequency Number of unobserved hH x hha hH x HH hH x hH locus A locus B allele (genotype) heterozygous sires heterozygous dams homozygous offspring Hhb hH hH HH hh hH HH 5:10278129 27:25728096 aB (aaBB) 105 2642 0.01 23 0 0 4382 4878 0 41 23 9:93731471 21:19601568 Ab (AAbb) 70 2296 0.009 22 0 0 4462 4626 0 34 22 11:32153221 21:19601568 ab (aabb) 68 1961 0.009 19 0 0 4155 4557 0 29 19 24:39650461 29:28523337 aB (aaBB) 80 2688 0.008 18 1 0 3337 3648 0 36 18 15:5766975 23:24441569 Ab (AAbb) 65 1650 0.007 18 0 0 3509 3619 0 28 18 15:54803685 16:79553324 Ab (AAbb) 59 1841 0.007 18 0 0 3583 3671 0 28 18 7:4436640 13:72570537 aB (aaBB) 72 2524 0.008 18 0 0 3598 3770 0 27 18 11:30575025 25:4667894 Ab (AAbb) 71 2506 0.009 18 0 0 3720 3792 0 29 18 8:82622674 20:65064787 Ab (AAbb) 102 1415 0.011 18 0 0 5954 5991 0 32 18 11 4:84118706 10:29375348 AB (AABB) 71 2744 0.009 17 0 0 3579 3909 0 31 17 1:86931672 16:12741497 ab (aabb) 69 1183 0.006 17 0 0 3835 3965 0 27 17 1:7782816 9:81277154 ab (aabb) 94 2809 0.009 17 0 0 4033 4274 0 30 17 7:26984291 12:18804912 Ab (AAbb) 70 1535 0.009 16 0 1 2638 2839 0 24 16 14:82345678 20:33773531 ab (aabb) 65 2504 0.01 16 0 0 2630 2874 0 24 16 6:47791866 25:34622211 Ab (AAbb) 66 2548 0.007 16 0 0 2925 3044 0 26 16 14:20606100 17:19346381 ab (aabb) 62 1998 0.007 16 0 0 3168 3284 0 25 16 9:7602607 10:57131487 aB (aaBB) 67 1900 0.007 16 0 0 3207 3399 0 24 16 10:44445350 17:19346381 aB (aaBB) 65 2085 0.007 16 0 0 3297 3399 0 25 16 10.3389/fgene.2023.1132796 15:62512689 16:79449472 AB (AABB) 58 1710 0.007 16 0 0 3511 3597 0 24 16 3:30194236 25:20619710 Ab (AAbb) 39 1595 0.006 16 0 0 3497 3630 0 27 16 a The carrying status of the artificial haplotype in the parents’ generation. b The carrying status of the artificial haplotype in the offspring generation; h: artificial haplotype; H: alternative haplotypes. frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 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, acanovas@uoguelph.ca. FIGURE 1 The observed, under- and over-represented offspring across matings in two-locus genotype with (A) additive-by-additive and (B) Ethics statement dominance-by-dominance epistatic effects. Ethical review and approval was not required for the animal study because the study was exempt from ethical approval procedures. et al., 2017) and conception rate at first service (Galliou et al., 2020) were also found for BAT4:85759993 and BAT21:37978844. For the Author contributions particular pair “BAT7:10835967 and BAT8:101251865” we find the existence of QTL for the interval from first to last insemination SI-L, JC, and AC conceived the research and developed the (Zhang et al., 2019) and daughter pregnancy rate (Cole et al., 2011). models. SI-L developed the program for data analyses and The SNP BAT7:10835967 also interacted with BAT10:29761646 and performed the corresponding analyses. SI-L wrote the first draft where QTL were reported for calving ease, daughter pregnancy rate of the manuscript. SI-L, JC, AC, FS, and FM participated in the and stillbirth by Cole et al. (2011), calving ease by Müller et al. interpretation of the results and revision of the manuscript. AC and (2017), interval to first estrus after calving by Liu et al. (2017) and JC supervised the research. AC secured the funding for the project. conception rate at first service by Galliou et al. (2020). Moreover, the All authors read and approved the final manuscript. 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 Funding epistatic TRD findings, however, further research investigating more in-depth the genomic regions with epistatic TRD Authors thank the financial support in main part by Ontario concerning their biological and functional implications are needed. 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 5 Conclusion and Forage Science Cluster (FDE.13.17) funded by the Canadian Beef Cattle Check-Off, Beef Cattle Research Council (BCRC), Our results aimed to elucidate the prevalence and patterns of Alberta Beef Producers, Alberta Cattle Feeders’ Association, Beef epistatic transmission ratio distortion across the Holstein genome. Farmers of Ontario, La Fédération des Productuers de bovins du Different epistatic TRD patterns were observed. Using genotypic Québec, and Agriculture and Agri-Food Canada’s Canadian Frontiers in Genetics 12 frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 Agricultural Partnership. This study was also supported by NSERC Publisher’s note (Natural Sciences and Engineering Research Council). All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated Acknowledgments organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or The authors acknowledge the data provided by Lactanet claim that may be made by its manufacturer, is not guaranteed or (Guelph, Ontario, Canada). endorsed by the publisher. Conflict of interest Supplementary material The authors declare that the research was conducted in the The Supplementary Material for this article can be found online absence of any commercial or financial relationships that could be at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1132796/ construed as a potential conflict of interest. full#supplementary-material References Arends, D., Karst, S., Heise, S., Korkuc, P., Hesse, D., and Brockmann, G. A. (2022). Crow, J. F. (1999). Unmasking a cheating gene. Science 283, 1651–1652. doi:10.1126/ Transmission distortion and genetic incompatibilities between alleles in a science.283.5408.1651 multigenerational mouse advanced intercross line. GENETICS 220 (1), iyab192. Cutter, A. D. (2012). The polymorphic prelude to Bateson-Dobzhansky-Muller doi:10.1093/genetics/iyab192 incompatibilities. Trends Ecol. Evol. 27, 209–218. doi:10.1016/j.tree.2011.11.004 Ballinger, M. A., and Noor, M. A. F. (2018). Are lethal alleles too abundant in Dobzhansky, T. H. (1937). Genetics and the origin of species. New York, NY: Columbia humans? 2018. Trends Genet. 34, 87–89. University Press. Bateson, W. (1909). “Heredity and variation in modern lights,” in Darwin and Fickel, J., and Weyrich, A. (2011). “Female mate choice in rodents,” in From genes to modern science. Editor A. C. Seward (Cambridge, MA, United Kingdom: Cambridge animal behavior. Editors M. Inoue-Murayama, S. Kawamura, and A. Weiss (Japan: University Press), 85–10. Springer Japan), 3–33. doi:10.1007/978-4-431-53892-9_1 Behrouzi, P., and Wit, E. C. (2018). Detecting epistatic selection with partially Fishman, L., Aagaard, J., and Tuthill, J. C. (2008). Toward the evolutionary genomics observed genotype data by using copula graphical models. J. R. Stat. Soc. C 68, of gametophytic divergence: Patterns of transmission ratio distortion in monkeyflower 141–160. doi:10.1111/rssc.12287 (mimulus) hybrids reveal a complex genetic basis for conspecific pollen precedence. Beissinger, T. M., Gholami, M., Erbe, M., Weigend, S., Weigend, A., de Leon, N., et al. Evolution 62, 2958–2970. doi:10.1111/j.1558-5646.2008.00475.x (2015). Using the variability of linkage disequilibrium between subpopulations to infer François, O., Martins, H., Caye, K., and Schoville, S. D. (2016). Controlling false sweeps and epistatic selection in a diverse panel of chickens. Heredity 116, 158–166. discoveries in genome scans for selection. Mol. Ecol. 25 (2), 454–469. doi:10.1111/mec. doi:10.1038/hdy.2015.81 13513 Borowsky, R., and Cohen, D. (2013). Genomic consequences of ecological speciation Galliou, J. M., Kiser, J. N., Oliver, K. F., Seabury, C. M., Moraes, J. G. N., Burns, G. W., in Astyanax cavefish. PLoS ONE 8 (11), e79903. doi:10.1371/journal.pone.0079903 et al. (2020). Identification of loci and pathways associated with heifer conception rate in Boschiero, C., Moreira, G. C. M., Gheyas, A. A., Godoy, T. F., Gasparin, G., Mariani, P. U.S. Holsteins. U.S. holsteins. Genes. 11, 767. doi:10.3390/genes11070767 D. S. C., et al. (2018). Genome-wide characterization of genetic variants and putative Giesbers, A. K. J., den Boer, E., Ulen, J. J. W. E. H., van Kaauwen, M. P. W., Visser, R. regions under selection in meat and egg-type chicken lines. BMC Genomics 25 (1), 83. G. F., Niks, R. E., et al. (2019). Patterns of transmission ratio distortion in interspecific doi:10.1186/s12864-018-4444-0 Lettuce hybrids reveal a sex-Independent gametophytic barrier. Genetics 211, 263–276. Brennan, A. C., Hiscock, S. J., and Abbott, R. J. (2014). Interspecific crossing and doi:10.1534/genetics.118.301566 genetic mapping reveal intrinsic genomic incompatibility between two Senecio species Gilbert-Diamond, D., and Moore, J. H. (2011). Analysis of gene-gene interactions. that form a hybrid zone on Mount Etna, Sicily. Heredity 113, 195–204. doi:10.1038/hdy. Curr. Protoc. Hum. Genet. 1-14. 2014.14 Haddad, R., Meter, B., and Ross, J. A. (2018). The genetic architecture of intra-species Bulmer, M. G. (1971). The effect of selection on genetic variability. Am. Nat. 105, hybrid mito-nuclear epistasis. Front. Genet. 9, 481. doi:10.3389/fgene.2018.00481 201–211. doi:10.1086/282718 Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and Carlborg, O., and Haley, C. S. (2004). Epistasis: Too often neglected in complex trait their applications. Biometrika 57, 97–109. doi:10.1093/biomet/57.1.97 studies? Nat. Rev. Genet. 5, 618–625. doi:10.1038/nrg1407 Hemani, G., Shakhbazov, K., Westra, H. J., Esko, T., Henders, A. K., McRae, A. F., Casellas, J., Cañas-Álvarez, J. J., González-Rodríguez, A., Puig-Oliveras, A., Fina, M., et al. (2014). Detection and replication of epistasis influencing transcription in humans. Piedrafita, J., et al. (2017). Bayesian analysis of parent-specific transmission ratio distortion Nature 508, 249–253. doi:10.1038/nature13005 in seven Spanish beef cattle breeds. Anim. Genet. 48, 93–96. doi:10.1111/age.12509 Hollander, W. F. (1955). Epistasis and hypostasis. J. Hered. 46, 222–225. doi:10.1093/ Casellas, J., Gularte, R. J., Farber, C. R., Varona, L., Mehrabian, M., Schadt, E. E., et al. oxfordjournals.jhered.a106562 (2012). Genome scans for transmission ratio distortion regions in mice. Genetics 191, 247–259. doi:10.1534/genetics.111.135988 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. Casellas, J., Id-Lahoucine, S., and Cánovas, A. (2020). Discriminating between allele- pone.0087330 and genotype-specific transmission ratio distortion. Anim. Genet. 51 (6), 847–854. doi:10.1111/age.13007 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 Casellas, J., Manunza, A., Mercader, A., Quintanilla, R., and Amills, M. (2014). A with transmission ratio distortion across the whole genome: TRDscan v.1.0. J. Dairy Sci. flexible Bayesian model for testing for transmission ratio distortion. Genetics 198, 102, 3175–3188. doi:10.3168/jds.2018-15296 1357–1367. doi:10.1534/genetics.114.169607 Id-Lahoucine, S., Casellas, J., Fonseca, P. A. S., Suárez-Vega, A., Schenkel, F. S., and Cole, J. B., Wiggans, G. R., Ma, L., Sonstegard, T. S., LawlorJr, T. J., Crooker, B. A., Cánovas, A. (2022). Deviations from mendelian inheritance on bovine X-chromosome et al. (2011). Genome-wide association analysis of thirty one production, health, revealing recombination, sex-of-offspring effects and fertility-related candidate genes. reproduction and body conformation traits in contemporary U.S. Holstein cows. Genes 13, 2322. doi:10.3390/genes13122322 BMC Genomics 12, 408. doi:10.1186/1471-2164-12-408 Id-Lahoucine, S., and Casellas, J. (2017). Impact of incomplete pedigree data and Corbett-Detig, R., Zhou, J., Clark, A. G., Hartl, D. L., and Ayroles, J. F. (2013). Genetic independent culling level pre-selection on the genetic evaluation of livestock: A incompatibilities are widespread within species. Nature 504, 135–137. doi:10.1038/ simulation study on lamb growth. Livest. Sci. 198, 76–81. doi:10.1016/j.livsci.2017. nature12678 02.011 Frontiers in Genetics 13 frontiersin.org Id-Lahoucine et al. 10.3389/fgene.2023.1132796 Id-Lahoucine, S. (2020). Mendelian inheritance in the genomics and big data era: Niedzicka, M., Dudek, K., Fijarczyk, A., Zielinski, P., and Babik, W. (2017). Transmission ratio distortion phenomenon in cattle breeds. [dissertation]. Guelph, Linkage map of Lissotriton newts provides insight into the genetic basis of Canada: University of Guelph. Available at: https://hdl.handle.net/10214/18142. reproductive isolation. G3 Genes, Genomes, Genet. 7, 2115–2124. doi:10.1534/ g3.117.041178 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- Ohta, T. (1982). Linkage disequilibrium due to random genetic drift in finite 019-38689-2 subdivided populations. Proc. Natl. Acad. Sci. USA. 79, 1940–1944. doi:10.1073/ pnas.79.6.1940 Jeffreys, H. (1984). Theory of probability. Oxford: Clarendon Press. Phillips, P. C. (2008). Epistasis – The essential role of gene interactions in the structure Kass, R. E., and Raftery, A. E. (1995). Bayes factors. J. Am. Stat. Assoc. 90, 773–795. and evolution of genetic systems. Nat. Rev. Genet. 9, 855–867. doi:10.1038/nrg2452 doi:10.1080/01621459.1995.10476572 Phillips, P. C. (1998). The language of gene interaction. Genetics 149, 1167–1171. Kempthorne, O. (1969). An introduction to genetic statistics. Ames Iowa: The Iowa doi:10.1093/genetics/149.3.1167 University Press. Sargolzaei, M., Chesnais, J. P., and Schenkel, F. S. (2014). A new approach for efficient Kiser, J. N., Keuter, E. M., Seabury, C. M., Neupane, M., Moraes, J. G. N., Dalton, J., genotype imputation using information from relatives. BMC Genomics 15, 478. doi:10. et al. (2019). Validation of 46 loci associated with female fertility traits in cattle. BMC 1186/1471-2164-15-478 genomics 20, 576. doi:10.1186/s12864-019-5935-3 Spiegelhalter, D. J., Best, N. G., Carlin, B. P., and van der Linder, A. (2002). Bayesian measures Leppälä, J., Bokma, F., and Savolainen, O. (2013). Investigating incipient speciation in of model complexity and fit. J. Roy. Stat. Soc. B 64, 583–639. doi:10.1111/1467-9868.00353 Arabidopsis lyrata from patterns of transmission ratio distortion. Genetics 194, 697–708. doi:10.1534/genetics.113.152561 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. Liu, A., Wang, Y., Sahana, G., Zhang, Q., Liu, L., Lund, M. S., et al. (2017). Genome- J. Anim. Sci. 91, 1562–1569. doi:10.2527/jas.2012-5584 wide association studies for female fertility traits in Chinese and nordic holsteins. Sci. Rep. 7, 8487. doi:10.1038/s41598-017-09170-9 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. Mackay, T. F. C. (2014). Epistasis and quantitative traits: Using model Editors D. S. Huang, K. Han, and A. Hussain (Switzerland: Springer), 21–32. organisms to study gene-gene interactions. Nat. Rev. Genet. 15, 22–33. doi:10. 1038/nrg3627 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. Maltecca, C., Weigel, K. A., Khatib, H., Cowan, M., and Bagnato, A. (2008). Whole- doi:10.1186/s13040-017-0143-7 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, Thomasen, J. R., Guldbrandtsen, B., Sorensen, P., Thomsen, B., and Lund, M. S. 27–34. doi:10.1111/j.1365-2052.2008.01793.x (2008). Quantitative trait loci affecting calving traits in Danish Holstein cattle. J. Dairy Sci. 91, 2098–2105. doi:10.3168/jds.2007-0602 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. Turelli, M., and Orr, H. A. (2000). Dominance, epistasis and the genetics of 1093/genetics/143.4.1739 postzygotic isolation. Genetics 154, 1663–1679. doi:10.1093/genetics/154.4.1663 Muller, H. J. (1942). Isolating mechanisms, evolution and temperature. Biol. Symp. 6, van Boven, M., and Weissing, F. J. (2001). Competition at the mouse t complex: Rare alleles 71–125. are inherently favored. Theor. Popul. Biol. 60, 343–358. doi:10.1006/tpbi.2001.1551 Müller, M. P., Rothammer, S., Seichter, D., Russ, I., Hinrichs, D., Tetens, J., et al. VanRaden, P. M., Olson, K. M., Null, D. J., and Hutchison, J. L. (2011). Harmful (2017). Genome-wide mapping of 10 calving and fertility traits in Holstein dairy cattle recessive effects on fertility detected by absence of homozygous haplotypes. J. Dairy Sci. with special regard to chromosome 18. J. Dairy Sci. 100, 1987–2006. doi:10.3168/jds. 94, 6153–6161. doi:10.3168/jds.2011-4624 2016-11506 Zhang, Z., Kargo, M., Liu, A., Thomasen, J. R., Pan, Y., and Su, G. (2019). Genotype- Neel, J. V., and Schull, W. J. (1954). Human heredity. Chicago: University of Chicago by-environment interaction of fertility traits in Danish Holstein cattle using a single- Press. step genomic reaction norm model. Heredity 123, 202–214. Frontiers in Genetics 14 frontiersin.org