Although many studies attempted to obtain estimates of nonadditive genetic effects from the total genetic variance in cattle and pigs [4, 13–14, 17, 26–27], the data size used in these studies were relatively small, possibly hampering their identification. Additionally, the number of traits studied was limited. In this study, a larger number of pigs were used (26,735 to 27,647 pigs and 5,338 sows with both genotypes and phenotypes for the production and reproduction traits). In addition, in the current study additive-by-dominance and dominance-by-dominance effects were also investigated. A similar approach was taken by Vitezica et al. [12]. However, the data size used in their study was relatively small (3,619 sows), and only one trait was analyzed (litter size).
Variance decomposition with additive, dominance, and epistatic effects
In general, in the current study dominance genetic variance explained a small proportion of total variation (0 for reproductive traits and 1-2% for production traits), smaller than estimates in pigs reported by Guo et al. [4], Xiang et al. [27], and Norris et al. [28] (4%-10% for NBA and 5% for daily gain and BF). Duroc pigs were used by Guo et al. [4] and Norris et al. [28], while Landrace and Yorkshire pigs were used by Xiang et al. [27]. However, these studies did not include the epistatic effects. Same as our study, Vitezica et al. [12] included the additive-by-additive, additive-by-dominance, and dominance-by-dominance effects in the statistical model for liter size and observed about 1.9% of dominance genetic variance, which is much lower than the estimates reported by the above studies but still larger than the estimates by our research. However, in Vitezica et al. [12] the breed used for the study was not specified, making it challenging to compare results from the current research. Regardless, the above results suggest that excluding the epistatic effects from the statistical model could result in biased estimates of dominance genetic variance. Vitezica et al. [12] also compared a full model that included the epistatic effect versus the reduced model that did not include the epistatic effect using the same data. Indeed, they found that excluding the epistatic effects generated a larger estimate of dominance genetic variance. The same was true for the additive genetic effect, estimates of additive genetic variance were biased upward when excluding the dominance and epistatic effects from the model. This is in agreement with our study (Table 1). In the same paper, Vitezica et al. [12] discussed how variance component estimates could differ if the effects in the model are not orthogonal. The assumption of orthogonality depends on HWE and does not hold in our study. This might have affected our estimates to some degree. This could also explain why the additive effect variance increased when nonadditive effects were excluded from the model, especially for NBA. Joshi et al. [29] used both the method that assumed HWE and the natural and orthogonal interaction approach that did not require HWE for constructing genomic relationship matrices. They found that estimates of variance components were only marginally different whether HWE was attained or not. Still, this had little to no impact on the qualitative outcomes related to the genetic architecture of the trait. The same was true for the traits investigated in our study.
In our study, the proportion of variation explained by epistatic effects was large for reproductive traits, in contrast to what was found by Vitezica et al. [12]. A large percentage of additive-by-additive effect variances were found for NBA and NW (31% and 12%, respectively), while the dominance-by-dominance variance was the largest for NBD (22%). Estimates of the epistatic effects for reproductive traits in pigs are scarce in literature. Vitezica et al. [12] reported a very low percentage of epistatic effect (4%) for littler size. However, no specific information on the pig breed was provided in their study. To further validate our results, two additive-by-additive genomic relationship matrices were compared using two approaches proposed by Jiang and Reif [9] and Vitezica et al. [11], respectively. The method proposed by Vitezica et al. [11] included interactions of loci with themselves (intra-locus interaction), which was excluded by Jiang and Reif [9]. Based on our data, the models using two different additive-by-additive genomic relationship matrices resulted in very similar estimates of variance components (Table 2). This suggests that when the number of markers is large enough, the effect of intra-locus interactions can be ignored. It is normally accepted that fitness-related traits display relatively high nonadditive variance [30]. Litter size is an important aspect of fitness in mammals and is lowly heritable in pigs [31–33]. Nonadditive genetic effects could be partially responsible for low heritability estimates observed in combination with high total genetic variance. Our studies found a large proportion of nonadditive genetic variance, particularly the epistatic genetic variances for NBA and NW. Literature in swine is scarce, however, Peripato et al. [34] observed that QTL and their interactions explain nearly 49% of the phenotypic variation for litter size in mice.
Table 2
Estimates of inbreeding depression and proportions of variance components in Duroc pigs using alternative genomic additive-by-additive relationship matrices
Trait
|
Sample size
|
b (P-value)
|
A (P-value)
|
D (P-value)
|
AA (P-value)
|
p (P-value)
|
NBA1
|
3246
|
-21.26 (< 0.001)
|
0.05 (0.01)
|
0 (1)
|
0.31 (< 0.001)
|
-
|
NBA2
|
3246
|
-21.26 (< 0.001)
|
0.05 (0.01)
|
0 (1)
|
0.31 (< 0.001)
|
-
|
NBD1
|
3246
|
-3.25 (0.67)
|
0.03 (0.05)
|
0 (1)
|
0.07 (0.36)
|
-
|
NBD2
|
3246
|
-3.25 (0.67)
|
0.03 (0.05)
|
0 (1)
|
0.07 (0.36)
|
-
|
NW1
|
3246
|
-14.06 (0.07)
|
0.04 (0.02)
|
0 (1)
|
0.10 (0.18)
|
-
|
NW2
|
3246
|
-14.06 (0.07)
|
0.04 (0.02)
|
0 (1)
|
0.10 (0.18)
|
-
|
BF1
|
26189
|
-15.03 (< 0.001)
|
0.34 (< 0.001)
|
0.01 (< 0.001)
|
0.11 (< 0.001)
|
0.06 (< 0.001)
|
BF2
|
26189
|
-15.03 (< 0.001)
|
0.34 (< 0.001)
|
0.01 (< 0.001)
|
0.11 (< 0.001)
|
0.06 (< 0.001)
|
WT1
|
26813
|
-22.81 (< 0.001)
|
0.08 (< 0.001)
|
0.02 (< 0.001)
|
0.16 (< 0.001)
|
0.11 (< 0.001)
|
WT2
|
26813
|
-22.81 (< 0.001)
|
0.08 (< 0.001)
|
0.02 (< 0.001)
|
0.16 (< 0.001)
|
0.11 (< 0.001)
|
MS1
|
26189
|
-11.26 (< 0.001)
|
0.18 (< 0.001)
|
0.01 (< 0.001)
|
0.11 (< 0.001)
|
0.05 (< 0.001)
|
MS2
|
26189
|
-11.26 (< 0.001)
|
0.18 (< 0.001)
|
0.01 (< 0.001)
|
0.11 (< 0.001)
|
0.05 (< 0.001)
|
1 used the additive-by-additive matrix proposed by Vitezica et al. (2017), 2 used the additive-by-additive matrix proposed by Jiang and Reif (2015). b: inbreeding depression; A: additive; D: dominance; AA: additive-by-additive; p: common litter environmental effect; NBA: number born alive; NBD: number born dead; NW: number weaned; BF: off-test back fat; WT: off-test weight; MS: off-test loin muscle depth |
Interestingly, large epistatic genetic variances were also observed for the production traits despite large additive effects being found, especially for BF and MS. Su et al. [8] reported that the additive-by-additive genetic effect accounted for 9.5% of total phenotypic variance for growth rate in Duroc pigs, mostly without genotype data. This estimate is very similar to WT (off-test weight, 10%) in our study, a trait related to growth rate. However, the additive-by-dominance effect, which was not included in Su et al. [8], having genetic variance much larger for WT (23%) than the additive-by-additive effect. Joshi et al. [29] observed a sizeable additive-by-additive genetic variance (17% of total phenotypic variance) for body weight at harvest in tilapia. Carcass trait BF showed large estimates of additive-by-additive (9%) and additive-by-dominance variance (9%), and MS showed large estimates of additive-by-dominance variance (17%) and dominance-by-dominance variance (10%). To our knowledge, no previous studies have reported estimates of epistatic effects on carcass traits in pigs. However, in beef cattle, Bolormaa et al. [17] found very low estimates of additive-by-additive variance for several carcass traits, which disagrees with our studies, in which 6% and 4% of additive-by-additive effects were observed for BF and MS in pigs, respectively. This may be because different species and slightly different carcass traits were used. Our findings also provide empirical evidence that nonadditive effects can be considerable for traits with relatively higher heritability, e.g., carcass traits.
Inbreeding depression is widely assumed to be deleterious for economically important traits, thus, can have a significant impact on animal production [35–37]. In the present study, NBA and more prevalently production traits showed a substantial inbreeding depression (Table 1), consistent with results reported by Christensen et al. [35] in pigs. Additionally, Xiang et al. [27] proved analytically that inclusion of genomic inbreeding as a covariate is necessary to obtain correct estimates of dominance variance in the presence of directional dominance. Vitezica et al. [12] also confirmed that genomic models that include nonadditive effects must also consider inbreeding depression to avoid upward bias of estimates of dominance variance for litter size in pigs. This agrees with our results. As shown in Table 3, a consistent reduction of estimates of dominance variance was observed for NBA and production traits, especially for WT, when inbreeding coefficient was included as a fixed covariate in the model. However, the effect of fitting genomic inbreeding on estimates of epistatic variance components for production traits in pigs was previously unknown. For the first time, our results show that the epistatic effects for production traits are largely unaffected by the inclusion of genomic inbreeding. While slight changes of estimates are observed these are likely due to random errors.
Table 3
Estimates of proportions of variance components in Duroc pigs with and without inbreeding coefficient
Trait
|
Sample size
|
b (P-value)
|
A (P-value)
|
D (P-value)
|
AA (P-value)
|
AD (P-value)
|
DD (P-value)
|
p (P-value)
|
NBA (W)
|
3246
|
-21.26 (0.01)
|
0.05 (0.01)
|
0 (1)
|
0.31 (< 0.001)
|
-
|
-
|
-
|
NBA (WO)
|
3246
|
-
|
0.05 (0.01)
|
0.004 (0.83)
|
0.30 (< 0.001)
|
-
|
-
|
-
|
NBD (W)
|
3246
|
-3.48 (0.75)
|
0.03 (0.14)
|
0 (1)
|
0.06 (0.50)
|
0 (1)
|
0.22 (0.55)
|
-
|
NBD (WO)
|
3246
|
-
|
0.03 (0.14)
|
0 (1)
|
0.06 (0.48)
|
0 (1)
|
0.20 (0.57)
|
-
|
NW (W)
|
3246
|
-13.95 (0.11)
|
0.03 (0.10)
|
0 (1)
|
0.12 (0.2)
|
0 (1)
|
-
|
-
|
NW (WO)
|
3246
|
-
|
0.03 (0.10)
|
0 (1)
|
0.12 (0.19)
|
0 (1)
|
-
|
-
|
BF (W)
|
26189
|
-15.35 (< 0.001)
|
0.33 (< 0.001)
|
0.008 (0.001)
|
0.09 (< 0.001)
|
0.09 (0.003)
|
0 (1)
|
0.05 (< 0.001)
|
BF (WO)
|
26189
|
-
|
0.33 (< 0.001)
|
0.011 (< 0.001)
|
0.09 (< 0.001)
|
0.08 (< 0.001)
|
0 (1)
|
0.05 (< 0.001)
|
WT (W)
|
26813
|
-23.71 (< 0.001)
|
0.08 (< 0.001)
|
0.02 (< 0.001)
|
0.10 (< 0.001)
|
0.23 (< 0.001)
|
0.02 (0.64)
|
0.10 (< 0.001)
|
WT (WO)
|
26813
|
-
|
0.08 (< 0.001)
|
0.03 (< 0.001)
|
0.10 (< 0.001)
|
0.22 (< 0.001)
|
0.02 (0.64)
|
0.10 (< 0.001)
|
MS (W)
|
26189
|
-12.27 (< 0.001)
|
0.19 (< 0.001)
|
0.007 (0.01)
|
0.05 (< 0.001)
|
0.17 (< 0.001)
|
0.10 (0.02)
|
0.04 (< 0.001)
|
MS (WO)
|
26189
|
-
|
0.19 (< 0.001)
|
0.01 (< 0.001)
|
0.06 (< 0.001)
|
0.17 (< 0.001)
|
0.10 (0.02)
|
0.04 (< 0.001)
|
b: inbreeding depression; A: additive; D: dominance; AA: additive-by-additive; AD: additive-by-dominance; DD: dominance-by-dominance; p: common litter environmental effect; W: statistical model includes inbreeding coefficient; WO: statistical model does not include inbreeding coefficient; NBA: number born alive; NBD: number born dead; NW: number weaned; BF: off-test back fat; WT: off-test weight; MS: off-test loin muscle depth |
Prediction performance
Epistasis did not improve the accuracy of prediction of litter size phenotype in pigs in Vitezica et al. [12]. In contrast, our results show that the addition of dominance and epistatic effects significantly improve the genomic prediction accuracy for NBA (Figure 1). Notably, a large percentage of additive-by-additive effect variance was observed in our study, compared to the 1.6% observed by Vitezica et al. [12]. This may occur due to two reasons. First, different pig breeds were used (reproductive traits haven’t been under selection in Duroc pigs); second, the trait studied were different with total number of pigs born per litter used by Vitezica et al. [12], while our study analyzed total number of pigs born alive per litter. It is noteworthy that total number of pigs born may also include total number of pigs born dead (NBD), which had a much lower estimate of additive-by-additive effect variance than NBA (Table 1), as shown in our study. There were very few studies reporting the genomic prediction of production traits by including epistatic effects in pigs. Su et al. [8] reported that the inclusion of dominance and additive-by-additive effect in a genomic prediction model increased prediction accuracy for daily gain in Duroc pigs. This is consistent with our findings. Genomic prediction accuracy was significantly improved when dominance and epistatic effects were included in the prediction model for off-test weight (WT, Figure 1), a trait similar to daily weight gain. However, in our study, the additive-by-dominance effect also contributed to increased accuracy. To the best of our knowledge, no studies have reported the genomic prediction accuracy, including dominance and epistasis, for carcass traits in pigs. Our results suggest that genomic prediction accuracy for loin muscle depth (MS) can be significantly improved when dominance and epistasis are included in the prediction model (Figure 1). The improvement can be mostly attributed to large additive-by-dominance and dominance-by-dominance effects. Guo et al. [4] reported improved genomic prediction accuracy for backfat and daily gain in Danish Duroc pigs when dominance effect was included in the prediction model. However, epistatic effects were not included, which may have resulted in upward biases for dominance variance estimates. It is noteworthy that in our study, predictions using only additive effect or the sum of additive and dominance effect estimated from the full model did not improve the accuracy for any of the traits. Our results suggest that the significant improvement of prediction accuracy for NBA, WT, and MS is to be attributed mainly to the epistatic effects, however further research is needed to confirm our findings.
Application of genomic selection with epistatic effects
Inclusion of nonadditive effects in the genomic prediction model can increase the prediction accuracy of phenotypes. It can also be used to predict the performance of future mating between two individuals. Studies have shown that inclusion of dominance effects in the matting program for dairy cattle can increase rates of genetic gain [13–14]. The inclusion of dominance effect requires probabilities of the genotypes, e.g. A1A1, A1A2, and A2A2, for the combination of the two individuals and the estimates of the additive and dominance effects for the same marker [38]. However, the use of models that include more complex interactions such as epistasis in the matting program is more challenging because the predicted performance of a mate should integrate the predictive performance over all possible future genotypic configurations of the expected progeny. One solution to this is to select the top markers with significant epistatic interaction and incorporate these markers in the matting program. Carlborg et al. [39] and Carter et al. [40] showed that epistatic interactions between four loci mediated a considerably higher response to selection of growth in chicken than predicted by a single-locus model. Furthermore, Martini et al. [41] also showed that selecting just a subset of the largest epistatic interaction effects has the potential to improve predictive ability. However, further studies are needed to evaluate the epistatic effects between markers across the genome using our data.
Genome-wide association study of dominance effect
Several GWAS studies have been performed to identify the dominance QTL for reproduction and production traits in pigs [18, 42–43]. Specifically, Coster et al. [42] detected significant markers with dominance effect for NBA and litter size in pigs of unspecified breed. Lopes et al. [43] identified one significant marker with dominance effect for teat number in Landrace pigs. Yang et al. [18] further detected two significant markers associated with dominance effect on growth rate in large white pigs. However, these studies were based on relatively small data size, which may result in lower statistical power to detect the dominance signals. A larger phenotypic and genotypic dataset was used in our study. Additionally, both one-step and two-step GWAS strategies were applied in this study to validate the results (Figures 2 and 3). Our study suggested that the two-step GWAS approach was comparable to the one-step full mixed model. The two-step method used in this study was proved to be an efficient alternative to detect dominance effects when fast implementations of full mixed models are not available. Interestingly, a strong dominance signal was identified on SSC6 (135-136 Mb) for BF with high LD between SNPs in this QTL region (r2 = 0.42-0.75). The same trait was investigated by Yang et al. [18], who did not find any dominance signals associated with BF. This may be because, first, a much smaller dataset was used by Yang et al. [18] than by our study; second, different pig breeds were used, Yang et al. [18] studied the large white pigs, whereas Duroc pigs were used in our research. It is known that the large white pigs have been selected mainly for reproductive traits. In contrast, the Duroc breed has been selected for production traits for a long time, resulting in different genetic architecture between the two populations.
A candidate gene AK5 was harbored in the QTL region (135-136 Mb on SSC6) identified for BF. AK5 was found to be strongly associated with increased body fat and may be helpful in the treatment of obesity [44]. From our results, AK5 may also play a role in determining the back fat of Duroc pigs. Larger estimate of dominance effect than additive effect for the most significant maker for BF suggested overdominance for this gene. Furthermore, the QTL region identified for BF overlapped with previously identified QTL in pigs associated with back fat at the 10th and last rib, average backfat thickness, marbling score, and intramuscular fat content, indicating a significant association between this QTL and backfat.