A bi-dimensional genome scan for prolificacy traits in pigs shows the existence of multiple epistatic QTL

Background Prolificacy is the most important trait influencing the reproductive efficiency of pig production systems. The low heritability and sex-limited expression of prolificacy have hindered to some extent the improvement of this trait through artificial selection. Moreover, the relative contributions of additive, dominant and epistatic QTL to the genetic variance of pig prolificacy remain to be defined. In this work, we have undertaken this issue by performing one-dimensional and bi-dimensional genome scans for number of piglets born alive (NBA) and total number of piglets born (TNB) in a three generation Iberian by Meishan F2 intercross. Results The one-dimensional genome scan for NBA and TNB revealed the existence of two genome-wide highly significant QTL located on SSC13 (P < 0.001) and SSC17 (P < 0.01) with effects on both traits. This relative paucity of significant results contrasted very strongly with the wide array of highly significant epistatic QTL that emerged in the bi-dimensional genome-wide scan analysis. As much as 18 epistatic QTL were found for NBA (four at P < 0.01 and five at P < 0.05) and TNB (three at P < 0.01 and six at P < 0.05), respectively. These epistatic QTL were distributed in multiple genomic regions, which covered 13 of the 18 pig autosomes, and they had small individual effects that ranged between 3 to 4% of the phenotypic variance. Different patterns of interactions (a × a, a × d, d × a and d × d) were found amongst the epistatic QTL pairs identified in the current work. Conclusions The complex inheritance of prolificacy traits in pigs has been evidenced by identifying multiple additive (SSC13 and SSC17), dominant and epistatic QTL in an Iberian × Meishan F2 intercross. Our results demonstrate that a significant fraction of the phenotypic variance of swine prolificacy traits can be attributed to first-order gene-by-gene interactions emphasizing that the phenotypic effects of alleles might be strongly modulated by the genetic background where they segregate.


Background
In the last few years, a major breakthrough in the understanding of the genetic factors that shape complex traits has been the demonstration that, in several species, a nonnegligible fraction of the genetic variance is explained by epistatic interactions. The recent identification of multiple epistatic QTL controlling complex traits in mice [1][2][3][4], chickens [5,6], and in model organisms such as yeast [7] and Drosophila melanogaster [8,9] has been a major achievement in the understanding of the genetic nature of complex traits. In addition, the discovery that gene expression is modulated, amongst others, by a plethora of regulatory RNAs with diverse functions and properties has added a new and thick layer of complexity in the subsequent identification of the polymorphisms involved in these interactions, since many of them might reside in non-coding regions [10].
In domestic species, traits relying on reproductive physiology, such as prolificacy and fecundity, have a notable impact on the financial outcome of farming enterprises. In pigs, prolificacy is a complex trait that displays a low heritability and strong heterosis [11]. One-dimensional studies have reported the existence of several QTL affecting litter size in pigs [12][13][14][15][16]. However, only one of the reported QTL was significant on a genome-wide level (p < 0.05) [16], and there was a general lack of positional concordance amongst different genome scans [17]. More importantly, these QTL studies exclusively dissected the additive and dominance components of litter size, thus neglecting the analysis of epistatic interactions that, paradoxically, are expected to explain a substantial portion of genetic variation of reproductive traits [18]. In consequence, many unsolved questions concerning to the genetic architecture of pig prolificacy still remain to be answered. Which are the specific contributions of dominance and epistasis in modelling the phenotypic expression of this complex trait? If epistasis is important, which are the dimensions, geometry and intricacy of the network of interacting loci and which types of epistatic interactions are more relevant? In a cross between two inbred mice strains Peripato et al. [3] demonstrated the existence of eight interacting QTL that explain almost 49% of the phenotypic variance of litter size in this cross. These results highlighted the importance of non additive genetic variance as a fundamental component of prolificacy. Nevertheless, laboratory mice strains are usually bred in a very stable environment, where fluctuations are kept to a minimum, and they have been the subject of an intense process of genetic selection without parallel in any other mammal species. Moreover, mice belong to a different superorder (Euarchontoglires) than most of mammalian domestic species (Laurasiatheria), so it is reasonable to expect that in these two distantly related taxonomic groups the biology of reproduction can differ in many instances.
The relevance of the aforementioned questions led us to analyse the genetic architecture of prolificacy traits in pigs. In this way, we have performed an F 2 intercross between two distinct European and Asian breeds, the Iberian and Meishan porcine breeds. Chinese Meishan is one of the most prolific pig breeds of the world being an excellent candidate population to perform these kinds of studies [19]. Iberian is an autochthonous Spanish breed with a very low prolificacy [20]. There is a very marked phenotypic difference for prolificacy traits between these two breeds (around 7 piglets per parity), being 14.3 the mean for the number of piglets born alive per parity of the Meishan breed [19] and 7.0 the mean for this trait of the Iberian Guadyerbas strain [20]. Interestingly, the ancestors of these breeds are assumed to have diverged at least 150,000 years before present without subsequent introgressions [21]. In consequence, it is reasonable to expect that these breeds have evolved, since then, by following independent processes of artificial selection and genetic drift, thereby establishing different adaptive epistatic genetic complexes [22]. In the current work, we have performed both a onedimensional and a bi-dimensional genome-wide scans for prolificacy traits by employing this Iberian by Meishan F 2 intercross as a genetic resource. Our main objective was to elucidate if epistasis makes a major contribution in shaping the phenotypic variability of prolificacy in pigs.

Phenotypic data recorded in the F 2 sows and linkage map
A description of the data and statistics of phenotypic records of number of piglets born alive (NBA) and total number of piglets born (TNB) in the F 2 population is given in Table 1. The phenotypic variance was 10.24 and 9.61 for TNB and NBA, respectively. The linkage map of the 115 markers used in the QTL analyses is shown in Table 2. Marker order and distances as well as average chromosome lengths were in general agreement with other mapping projects and the USDA genome database http://www.animalgenome.org/pig/. Markers provided coverage of the 18 autosomes, with intervals between adjacent markers that were below 20 cM whenever possible. The average marker interval was 17.6 cM (sex-averaged map distance).

One dimensional genome scan for TNB and NBA
Results of the whole-genome scan using a single-QTL model for TNB and NBA are summarized in Table 3. Two genome-wide highly significant QTL were identified on SSC13 (P < 0.001) and SSC17 (P < 0.01) at similar positions for both traits. In SSC13, the single QTL for NBA and TNB were found at positions 50 and 55 cM, respectively (Table 3), sharing an overlapping region located between markers SW398 and SW2440. The most likely position for the QTL found on SSC17 was at 22 cM for both traits (Table 3).   Iberian allele was the one associated with an increase in 0.73 (± 0.19) piglets per copy for the QTL on SSC17. Moreover, this QTL on SSC17 displayed a negative dominance effect (-0.82 ± 0.29). We estimated the degree of dominance as the ratio d/a between the estimated dominance (d) and the absolute value of the additive effect (a).
Values of d/a larger than unity corresponds to overdominance, while a d/a ratio between 0 and 1 represents partial dominance. In both cases (QTL on SSC13 and SSC17), the estimated d/a values were consistent with a complete dominance situation. Similar values of additive and dominance effects were obtained for the TNB QTL on SSC13 and SSC17 (Table 3), a result that is not surprising since these two traits are highly correlated. Besides, the proportion of phenotypic variance explained by these two single QTL detected on SSC13 and SSC17 ranged from 2% to 3% for both TNB and NBA.

Identification of multiple interacting QTL for NBA and TNB
Results from the bi-dimensional genome scan for NBA are shown in Table 4. Four bi-dimensional genome-wide highly significant (P < 0.01) and five significant (P < 0.05) epistatic interactions between QTL were found. We confirmed that all the observed epistatic interactions were consistently detected across families rather than being the consequence of a single sire effect, a feature that is particularly important when the number of founder males is moderate or even small. The results obtained through the likelihood ratio test were further confirmed by using other approaches. First, a false discovery rate (FDR) was calculated based on the nominal P-values every 30 cM, with a result of 0.018 for a P-value lower than 0.001. Moreover, parametric bootstrapping confirmed the significance of the results.
A graphical overview of the epistatic interactions for NBA (red arrows) is shown in Figure 1. As much as twelve of the 18 pig autosomes (1,5,6,7,8,9,10,12,13,14,15 and 18) were involved in these interactions, forming a complex network with a non-radial geometry. This means that a specific region did not interact simultaneously with multiple loci, but with a very limited number of them (usually interactions were one to one). For example, the SSC12 region located at 11 cM, interacted significantly with another SSC12 region at 89 cM ( Figure 1; Table 4). Similarly, two non-overlapping SSC6 QTL regions showed epistatic interactions, one of them with QTL on SSC1 and SSC7 (SSC6, 54-69 cM) and another one with SSC14 (SSC6, 1-10 cM). As shown in Figure 1, other pig chromosomes exhibiting more than one significant interacting QTL were SSC1 (at positions 76 and 153 cM respectively), and SSC7 (at positions 28 cM, and 107 cM). An interesting feature of our analysis was that the highly significant NBA QTL identified in the one-dimensional genome scan (SSC13 at 50 cM and SSC17 at 22 cM) did not show any significant epistatic interaction with other regions across the genome, meaning that its mode of action is purely additive. In contrast, a NBA QTL found on another region of SSC13 (73 cM) had significant epistatic interactions with a QTL located at position 4 cM of SSC9. Similarly, in mice, Peripato et al. [3] identified two significant QTL for litter size in a one-dimensional genome scan (chromosomes 7 and 12) that did not emerge in the bidimensional analysis (chromosomes 2, 4, 5, 11, 14, 15 and 18). In the light of these results and ours, we could conclude that there is a low concordance between the QTL identified in one-and bi-dimensional genome scans. This means that, in general, the additive and epistatic components of prolificacy traits encompass different sets of genes.
With regard to the bi-dimensional genome-wide scan for TNB, we found three genome-wide highly significant (P < 0.01) and six significant (P < 0.05) epistatic interactions ( Table 5; Figure 1). Thirteen out of the 18 pig autosomes were involved in the epistatic QTL interactions for both traits. However, the network of interacting QTL for TNB was not identical to the one reported for NBA. In this sense, only around one third of the QTL epistatic interactions detected in the current study were coincident in both traits. This was an unexpected result since NBA and TNB display a high genetic correlation (r g~ 0.9) and they are expected to share a similar genetic architecture [11]. Moreover, the number of chromosomes displaying multiple interactions was higher for TNB (SSC1, 7, 12, 13 and 14) than for NBA. Notably, four epistatic QTL on SSC7 showed significant interactions with QTL located on SSC6 (60 cM), SSC13 (77 cM) and SSC1 (79 cM and 139 cM, respectively). Other epistatic interactions involved two Network representation of the epistatic QTL interactions in thirteen pig chromosomes (SSC) for prolificacy traits NBA (red arrows) and TNB (black arrows) Figure 1 Network representation of the epistatic QTL interactions in thirteen pig chromosomes (SSC) for prolificacy traits NBA (red arrows) and TNB (black arrows). 10 overlapping regions of SSC13 with one region of SSC7 and one SSC9 QTL, a feature that demonstrates the remarkable complexity and intricacy of these networks. Finally, it is worth mentioning that if we would have assumed a type 1 error α = 0.10, which corresponds to a genome-wide critical value of 25.14 for the LR test, we would have been able to detect 11 and 9 additional epistatic interactions between QTL across the genome for NBA and TNB, respectively (results not shown).

Partition of the phenotypic variance explained by epistatic QTL
We estimated the proportion of the phenotypic variance explained by epistatic genetic components for each TNB and NBA significant epistatic QTL pair. This contribution to the total phenotypic variance ranged from 3.26% (SSC14-SSC15) to 4.04% (SSC12-SSC12) for TNB and from 3.10% (SSC1-SSC6) to 3.62% (SSC9-SSC13) for NBA. The relative contribution of epistasis to the total phenotypic variance was estimated by adding the estimates of the partial epistatic effects of each epistatic interaction. The total phenotypic variance explained by the joint genetic effects of all epistatic QTL pairs for NBA and TNB was 37.6% and 42.4%, respectively. Nevertheless, we would like to mention that this approach might lead to an overestimation of the epistatic contribution because the models employed in our analysis are exclusive [23]. We have calculated the repeatability with a standard animal model, resulting in an estimated value of 0.27 for both TNB and NBA. This value should be considered as the limit of the total variance explained by genetic effects. This result evidences that our estimates of epistatic contributions are clearly overestimated but, in spite of this drawback, we think it is reasonable to assume that a significant proportion of the phenotypic variance of prolificacy is explained by the joint genetic effects of epistatic QTL. All four forms of epistasis (a × a, a × d, d × a and d × d) were detected among the significant epistatic pairs. In total, 9 a × a, 7 a × d, 11 d × a, and 16 d × d significant interactions were detected (Tables 4 and 5). For the SSC1-SSC6 (NBA) and SSC8-SSC14 (TNB) epistatic pairs, only one type of epistasis was significant (d × d, and a × d, respectively); whereas, the remaining pairs presented two or more types of epistasis, leading to more complex patterns of interactions. We plotted the genotypic values expected for each genotypic class (Figure 2), thus identifying 'dominance-by-dominance' and 'coadaptive' patterns of epistasis as described by Carlborg and Haley [18] and Carlborg et al. [4]. The 'dominance-by-dominance' pattern of epistasis was found in several epistatic QTL pairs and corresponded to interactions were a significant d × d effect was found. In positive d × d interactions (Figure 2a), the genotypic value of double heterozygotes I 1 M 1 -I 2 M 2 , (being I = Iberian alleles and M = Meishan alleles for loci 1 and 2) was superior to that of simple heterozygotes (I 1 M 1 -I 2 I 2 , I 1 M 1 -M 2 M 2 , I 1 I 1 -I 2 M 2 , and M 1 M 1 -I 2 M 2 ); a pattern that was reversed in negative d × d interactions (Figure 2b). According to Carlborg and Haley [18], the 'coadaptive' type of interaction occurs when positive a × a interaction effect is significant. This form of epistasis leads to enhanced performance of parental double homozygotes (I 1 I 1 -I 2 I 2 and M 1 M 1 -M 2 M 2 ) in comparison to the hybrid double homozygotes (I 1 I 1 -M 2 M 2 or M 1 M 1 -I 2 I 2 ). As an example, we found this specific pattern of epistasis in the NBA epistatic QTL pair SSC10-SSC15 (Figure 2c).

Classification of epistatic effects for NBA and TNB
Additionally, we found another epistatic pattern where only significant a × d and/or d × a effects were present. These pairs of interacting QTL were included in a category that we named 'additive-by-dominance' epistasis because a × d and d × a interactions showed exactly the same pattern. In positive 'additive-by-dominance' interactions, one locus is overdominant, neutral, or underdominant, depending on the genotype at the second locus. Conversely, the second locus is additive, with the favoured allele being dependent on the genotype at the first locus.
In negative 'additive-by-dominance' interactions, the genotypic values of homozygotes at one locus (_._I 2 I 2 and _._M 2 M 2 ) deviate from what is expectable under an additive model only when the second locus is heterozygous (I 1 M 1 -I 2 I 2 and I 1 M 1 -M 2 M 2 , Figure 2d).
Finally, we considered as a separate category those QTL pairs that could not be classified into any of the preceding epistasis models. The main common feature of QTL pairs gathered in this group was that they all had a significant a × d or d × a epistatic term together with significant a × a and/or d × d terms, leading to 'complex' patterns of interactions ( Figure 2e). These complex patterns might be explained by the existence of higher order interactions involving more than two loci or simply because small sample size for certain genotypic classes might lead to inaccurate estimates of the epistatic effects.

Identification of genome-wide significant QTL with additive and dominant effects on pig prolificacy
In the single QTL analyses, we found two genome-wide highly significant QTL located on SSC13 and SSC17 affecting NBA and TNB. So far, and to the best of our knowledge, only one study has described a genome-wide significant QTL (P < 0.05) detected at 88 cM on SSC15 for NBA [16]. In the current work, the mode of gene action strongly differed amongst QTL. For instance, the additive effect detected for the SSC13 QTL stemmed from the Meishan breed, and determined an increase of 0.71 and 0.61 piglets for NBA and TNB, respectively. This result is in agreement with the phenotypic differences observed between the purebred parental lines. Conversely, beneficial alleles for the QTL found on SSC17 appeared to derive from the Iberian breed, with increases of 0.73 and 0.68 piglets for NBA and TNB, respectively. These cryptic Iberian QTL alleles, which increase prolificacy, provide a compelling example of the complexity of the genetic architecture of these traits. In addition, significant dominance effects were also detected for the two mentioned QTL; however, the direction of the effects was complete dominance for the QTL on SSC13, and recessivity for the QTL at SSC17.
Our results confirm the existence of several previously reported suggestive QTL for litter size in pigs. De Koning et al. [13]

Phenotypic variation of prolificacy traits is affected by a complex network of epistatic QTL
The statistical evidence suggesting that epistasis might be an important component of reproductive traits in mice [3] led us to perform a bi-dimensional genome scan for NBA and TNB. This analysis allowed us to demonstrate that phenotypic variation of these traits can be strongly influenced by a complex network of interacting loci. In this sense, this is the first study that shows genome-wide significant epistatic QTL affecting prolificacy in pigs. After implementing the highly stringent Bonferroni correction  for multiple testing, as much as seven QTL interactions remained significant at the 1% bi-dimensional genomewide level, whereas 11 were significant at the 5% genomewide level. Rather than identifying one or a few master regions interacting with multiple loci (which would have yielded a star-like geometry of interactions), we found that most chromosomal regions interacted with a single counterpart. However, almost all chromosomes in which significant epistatic QTL were found participated in more than one interaction. Five of the significant epistatic interactions detected had pleiotropic effects on both TNB and NBA traits, whereas five were only related to TNB and three to NBA. These surprising differences for two highly genetically correlated traits [27], might be interpreted in the light of the statistical values of the Likelihood ratio obtained in the contrast between models and the stringent bi-dimensional genome-wide threshold assumed. For instance, the chromosome pair SSC10-SSC15 reached an LR value of 29.41 (p < 0.05) for NBA and an LR value of 22.33 ("no significant") for TNB. Alternatively, there might be a biological meaning behind the specific differences observed in the geometry of the NBA and TNB networks, thus indicating that although similar metabolic and physiological pathways may be implicated in the regulation of both traits, other mechanisms may be operating independently.

Epistatic QTL for pig prolificacy display different types of interactions
This genetic dissection of the epistatic component of prolificacy in pigs was completed with an analysis of the types of interactions detected in the bi-dimensional genomewide scan. All the epistatic QTL had at least one significant type of interaction (see tables 4 and 5). In total, nine pairs showed additive by additive (a × a) epistasis, which in certain circumstances can have a 'nullification effect' because epistasis might cancel out the effects of individual loci at intermediate frequencies, making difficult to detect them in a conventional one-dimensional genome scan [28]. Furthermore, seven pairs showed additive by dominance (a × d) epistasis, eleven pairs showed dominance by additive (d × a) epistasis, and sixteen pairs showed dominance by dominance (d × d) epistasis. All these forms of epistasis contribute to heterosis [29]. It has been widely supported that heterosis plays an important role in the genetic architecture of reproductive traits [11,12].
Several patterns of interactions among genotypes were identified and classified according to Carlborg and Haley [18] and Carlborg et al. [4]. Eight pairs were classified as 'dominance-by-dominance' and two pairs as 'coadaptive' epistasis. The remaining QTL showed different patterns of epistasis from those previously described in the literature and were, therefore, grouped into two additional categories that we named 'additive-by-dominance' and 'com-plex' epistasis. 'dominance-by-dominance' epistasis leads to a deviated performance for the double heterozygotes compared to single heterozygotes. Among the d × d interactions found in our study, six had positive sign ( Figure  2a) and two had negative sign (Figure 2b). In positive d × d interactions, the genotypic value of the double homozygotes is superior to the simple heterozygotes. Conversely, this pattern is reversed when the sign is negative. Coadaptive epistasis occurs when double homozygotes from the same parental line show enhanced performance [18]. As mentioned above, we only found two QTL pairs which displayed this form of epistasis (Figure 2c). Coadaptive epistasis is fundamental to interpret post-zygotic reproductive isolation [30][31][32]. In each parental population, selection may have acted leading to fixation of different alleles at the relevant loci regulating prolificacy in a way that statistical epistasis is not apparent in either population. Second-generation hybrids (F 2 ) exhibit combination of alleles at different loci that were not present in any of the parental breeds, leading to the disruption of "coadapted" gene pools and the appearance of new phenotypes.
The third category included two pairs where only a × d and d × a effects were significant. We called this group 'additive-by-dominance' because the a × d and d × a interactions show the same patterns but with the roles of the loci reversed ( Figure 2d). Finally, we found QTL pairs which show 'complex' patterns of interactions characterised by having a significant a × d or d × a term together with significant a × a and/or d × d (Figure 2e). These interactions yielded complex patterns in which the genotypic value for a given genotype at one locus drastically changes depending on the genotype at the second locus.

Conclusions
In summary, the bi-dimensional genome scan of an Iberian × Meishan F 2 intercross has allowed to demonstrate that the genetic architecture of pig reproduction is mostly built as a complex network of interacting genes rather than being explained by the sum of the additive effects of a yet to be defined number of loci. Individual epistatic loci have moderate effects on the phenotypic variance of prolificacy and they are distributed in multiple chromosomal locations. Moreover, they display several types of interactions that sometimes cannot be easily ascribed to well defined models, thus suggesting the existence of additional interacting loci. In the next years, the fine mapping and identification of the causal mutations that explain the segregation of epistatic QTL in pigs will be a daunting but fascinating task that will likely unveil many of the secrets that underlie the biological grounds of complex traits. With regard to the two-QTL analyses, two different models were employed. The first model included the effects of two non interacting QTL. The statistical mixed model was:

Experimental design and phenotypic data
where a 1 and a 2 were the additive effects, d 1 and d 2 were the dominance effects for QTL 1 and 2, respectively. The coefficients c a1 , c d1 , c a2 and c d2 were calculated as before for locations 1 and 2.
The second model allowed for epistasis, i.e.: where The two-QTL analyses were performed using a full bidimensional genome scan. LR-tests comparing the models with and without the epistatic interaction effects (model 3 vs model 2) were computed at 1 cM intervals along the 2,017 cM of the 18 autosomes for each of the two QTL, leading to a total of 2,069,595 regression analyses for both NBA and TNB. The values of h 2 and p 2 used in this analysis were identical to those considered in model 1. The statistical contrast between models for evidence of epistasis was carried out using an LR test with 4 degrees of freedom in the numerator. As before, bi-dimensional genome-wide levels of significance were calculated using a Bonferroni correction assuming statistical independence every 30 cM. The genome-wide critical values of LR test for level of significance associated with type I errors α = 0. 10