Investigating the impact of non-additive genetic effects in the estimation of variance components and genomic predictions for heat tolerance and performance traits in crossbred and purebred pig populations

Abstract

Background

Non-additive genetic effects are often ignored in livestock genetic evaluations. However, fitting them in the models could improve the accuracy of genomic breeding values. Furthermore, non-additive genetic effects contribute to heterosis, which could be optimized through mating designs. Traits related to fitness and adaptation, such as heat tolerance, tend to be more influenced by non-additive genetic effects. In this context, the primary objectives of this study were to estimate variance components and assess the predictive performance of genomic prediction of breeding values based on alternative models and two independent datasets, including performance records from a purebred pig population and heat tolerance indicators recorded in crossbred lactating sows.

Results

Including non-additive genetic effects when modelling performance traits in purebred pigs had no effect on the residual variance estimates for most of the traits, but lower additive genetic variances were observed, especially when additive-by-additive epistasis was included in the models. Furthermore, including non-additive genetic effects did not improve the prediction accuracy of genomic breeding values, but there was animal re-ranking across the models. For the heat tolerance indicators recorded in a crossbred population, most traits had small non-additive genetic variance with large standard error estimates. Nevertheless, panting score and hair density presented substantial additive-by-additive epistatic variance. Panting score had an epistatic variance estimate of 0.1379, which accounted for 82.22% of the total genetic variance. For hair density, the epistatic variance estimates ranged from 0.1745 to 0.1845, which represent 64.95–69.59% of the total genetic variance.

Conclusions

Including non-additive genetic effects in the models did not improve the accuracy of genomic breeding values for performance traits in purebred pigs, but there was substantial re-ranking of selection candidates depending on the model fitted. Except for panting score and hair density, low non-additive genetic variance estimates were observed for heat tolerance indicators in crossbred pigs.

Background

Phenotypes are the result of the individual genotypic value, environmental deviation, and interactions between genotype and the environment. The genotypic values can be divided in additive (breeding values) and non-additive genetic effects due to interactions between alleles at the same locus (dominance) and two or more loci (epistasis) [1]. Non-additive genetic effects are often ignored in livestock genetic evaluations due to the greater importance of additive effects for breeding value estimation, and the fact that breeding animals pass alleles and not genotypes to their offspring [2,3,4,5]. Furthermore, non-additive genetic effects can be confounded with non-genetic effects, such as maternal and permanent environmental effects, and the analyses incorporating non-additive genetic effects are often more complex and computationally demanding [4,5,6].

Although breeders are mainly interested in the estimates of breeding values, fitting non-additive genetic effects in the genomic prediction models can increase the accuracy of the breeding values [7,8,9]. Furthermore, fitting non-additive genetic effects in the models can yield more accurate predictions of future performance and better mating plans to maximize offspring performance by exploring the combination of additive and non-additive genetic effects [2]. The increasing availability of large-scale genomic datasets has enabled the study of the non-additive actions of genes in numerous complex traits, including the evaluation of non-additive genetic effects in genomic prediction of breeding values (e.g., [2, 9, 10]), genome-wide association studies (e.g., [11,12,13]), and the development of new methods for the estimation and applications of non-additive genetic effects in breeding programs [6, 9, 14,15,16]. It is also important to note that accounting for inbreeding in variance component estimation and genomic prediction ensures unbiased variance components, minimizes bias issues in the estimates, accounts for directional dominance, and prevents inflation of dominance variance estimates [17,18,19].

Non-additive genetic effects are expected to have a greater impact in crossbred performance as compared to purebred populations, mainly due to heterosis [20,21,22]. As pig production is based on crossbreeding schemes, it is important to evaluate the magnitude of non-additive effects for different traits measured on purebred and crossbred populations, as well as evaluate the impact of these effects on the predictive ability of genomic breeding values. In general, adaptation and fitness traits tend to be more influenced by non-additive genetic effects (e.g., [17, 22,23,24]). In this context, heat tolerance indicators are key fitness traits to be included in swine breeding programs for improving pig health, welfare, and production [25, 26]. Breeding for improved heat tolerance is becoming more important due to the increase in global temperatures that directly impact swine health, performance, and welfare [27, 28]. Although estimates of additive genetic variance and heritability have been reported for heat stress traits in pigs (e.g., [29,30,31]), the impact of non-additive effects on direct indicators of heat tolerance in crossbred pigs is still unknown.

Non-additive genetic effects can be of different magnitude in purebred and crossbred populations [20,21,22], as well for across traits with different genetic architecture. Therefore, there is a need to evaluate the impact of non-additive genetic effects into genomic prediction models for a diverse set of traits, including production, fitness, and adaptation, within both purebred and crossbred populations. In this context, the primary objectives of this study were to estimate additive and non-additive genetic variance components for different trait groups (i.e., production, fitness, and adaptation traits) in independent purebred and crossbred populations. In addition, the impact of non-additive effects in the accuracy of genomic prediction of purebred animals was evaluated.

Results

The results obtained are presented separately for each of the two datasets evaluated.

Dataset 1: purebred pig population

Variance components and heritability estimates

Variance components and heritabilities were estimated for five traits representing phenotypes routinely collected from birth to breeding in a single nucleus pig line population [32]. Table 1 presents the estimates of variance components and Table 2 shows the estimates of heritability and non-additive genetic variance ratios for all five traits in the purebred pig population. All the fitted models converged for all traits, with the exception of the model including all non-additive genetic effects evaluated and inbreeding depression (MAIDE3) that did not converge for one of the traits (T3). When inbreeding was included as a covariate in the models, a small decrease in the additive genetic variance was observed (Table 1). For most of the traits, lower residual and additive genetic variance estimates were observed when incorporating non-additive genetic effects in the analyses.

Table 1 Variance component estimates based on models including or not inbreeding and non-additive genetic effects

Table 2 Heritability estimates based on models including or not inbreeding and non-additive genetic effects

Low dominance variance components were estimated, ranging from 0 to 0.0576 for T1, T2, T3, and T4, and ranging from 138.0280 to 171.8660 for T5. Notably, the dominance variance ratios across all traits were found to be small, accompanied by relatively large standard errors. T5 had the highest dominance variance ratio, which ranged from 0.0397 to 0.0437, while for the other traits it ranged from 0 to 0.0206. However, the estimates of additive-by-additive epistatic variance were notably higher than the dominance variance estimates. These ratios were approximately 0.12 for T1, 0.24 for T2, 0.28 to 0.43 for T3, 0.09 for T4, and 0.13 to 0.15 for T5, indicating a more substantial contribution of additive-by-additive epistasis to the variance components for these traits.

In general, there was a greater decrease in the additive genetic variance when additive-by-additive epistasis was included in the models (MAE, MAIE, MADE1, MAIDE1, MADE2, MAIDE2, MADE3, MAIDE3). However, for T1 the additive genetic variance estimates also presented larger standard errors (Table 1), indicating that the inclusion of epistasis in the model decreased precision of the additive genetic variance estimates for T1. As a consequence of the decrease in the additive genetic variance, a reduction in the heritability estimates (${h}_{a}^{2}$) were observed when additive-by-additive epistasis was fitted in the models (Table 2). Nevertheless, the variance components for additive-by-dominance and dominance-by-dominance epistasis were small and their variance ratios (${h}_{ad}^{2}$ and ${h}_{dd}^{2}$) were close to zero for almost all traits, except T3.

Predictive ability of genomic breeding values

Figure 1 presents the average bias, dispersion, and accuracy values obtained for the genomic breeding values for all five traits. As expected, greater accuracies were obtained for traits with higher ${h}_{a}^{2}$. Fitting inbreeding and/or dominance in the models did not impact the accuracy of breeding values compared to those obtained based on the MA model (Fig. 1). Models including epistasis (MAE, MAIE, MADE1, MAIDE1, MADE2, MAIDE2, MADE3, and MAIDE3) presented on average a relative decrease in accuracy of 31.58%, 13.97%, 3.50%, 0.82%, and 1.34% for T1, T2, T3, T4, and T5, respectively, when compared to MA (Fig. 1). The observed bias was close to zero and the dispersion estimates were close to one for all traits and models (Fig. 1), indicating that the inclusion of non-additive genetic effects in the model did not influence the dispersion of breeding values.

Fig. 1

Bias, dispersion, and accuracy estimates of genomic breeding values obtained using the Linear Regression method

In this section, all models incorporating non-additive genetic effects were compared to the additive animal model. To assess the impact of including non-additive genetic effects in the model on animals’ selection, we calculated the change in an animals’ ranking in comparison to the additive model (MA) based on Spearman’s rank correlation, as shown in Table 3. The percentage of commonly selected individuals when selecting the top 1%, 5%, 10%, and 20% based on different genomic prediction models were compared to the selected individuals based on MA. Table 4 shows the percentage of coincidence of the top 1%, 5%, 10%, and 20% best animals between the MA and the models including non-additive genetic effects. The models including epistasis effect presented the lowest percentage of coincidence of the top animals in comparison to the MA model. This percentage varied from 68.57% for trait T3 to 97.59% for T4. Although there was little change in the Spearman correlations across models (Table 3), the percentage of coincidence between the best animals according to MA were higher for models including inbreeding and/or dominance (MAI, MAD, and MAID). However, the percentage of coincidence decreased for models including epistasis (MAE, MAIE, MADE1, MAIDE1, MADE2, MAIDE2, MADE3, and MAIDE3), as shown in Table 4.

Table 3 Spearman’s rank correlation between the additive genetic model (MA) and other models considering non-additive genetic effects

Table 4 Percentage of commonly selected animals between the additive model and models including non-additive genetic effects

The values of Akaike Information Criterion (AIC) (Table S1) suggested that the inclusion of epistasis in the prediction models for T2, T3, and T5 improved the predictive performance of the models. Although, the Likelihood Ratio Test (LRT; Table S2) for the models fitted for the purebred pig population showed that the inclusion of epistasis was significant for the traits T2, T3 and T5, and incorporating dominance effect was also significant for T5.

Dataset 2: crossbred pig population

Variance components and heritability estimates

We also estimated variance components and heritabilities for heat tolerance indicators in a crossbred (Large White x Landrace) population. Table 5 presents variance component estimates and Table 6 presents the estimates of heritability and non-additive genetic variance ratios for traits related to heat stress response in lactating sows for models including or not non-additive genetic effects. All heat stress related traits have substantial additive genetic variability (Table 5) with heritability estimates ranging from 0.0231 to 0.2639 (Table 6). However, all traits had small non-additive genetic variance with large standard error estimates, except panting score (PS) and hair density (HD) (Table 5). PS presented a ${h}_{aa}^{2}$ estimate of 0.1342, which corresponds to 82.22% of the total genetic variance for both models including epistasis (MAIEpe and MAIDEpe). HD presented ${h}_{aa}^{2}$ estimates of 0.4232 and 0.4493 for MAIE and MAIDE models, respectively. The proportion of the total genetic variance explained by additive-by-additive epistasis for HD ranged from 64.95 to 69.59% for the MAIE and MAIDE models, respectively (Table 6).

Table 5 Variance components for traits related to heat tolerance based on models fitting non-additive genetic effects

Table 6 Heritability estimates and variance ratios considering models fitting non-additive genetic effects for heat tolerance indicators

As observed for the purebred population (Dataset 1), the inclusion of non-additive genetic effects in the model did not affect the residual variance estimates for most of the traits, except for HD (Table 5). However, for HD and PS there was also a decrease in the additive genetic variance estimates (Table 5) when epistasis was included in the model. When both dominance and epistasis were included in the repeatability model (for all heat stress-related traits, except HD), there was a decrease in the permanent environmental variance estimate for most of the traits.

Model comparison for the crossbred pig dataset

For the crossbred population, the models also presented similar values of Akaike Information Criterion (AIC) (Table S3). LRT (Table S4) for the models fitted for the crossbred pig population showed that fitting epistasis in the models had a significant effect only for PS. The models for the other traits did not have better goodness of fit when dominance and/or epistasis were fitted.

Discussion

We aimed to estimate variance components and evaluate the predictive ability of genomic models including non-additive genetic effects for different traits in purebred and crossbred populations. In general, we observed small estimates for dominance variance and the inclusion of additive-by-additive epistasis in the model reduced the estimates of additive genetic variance. The inclusion of non-additive genetic effects in models for genomic prediction resulted in changes in animals’ ranking and in the groups of animals to be selected. Since two different independent populations were used in this study, the Discussion section was structured into two subsections to discuss the findings for each population separately.

Dataset 1: purebred pig population

One way to account for inbreeding depression in genomic prediction is by fitting the inbreeding coefficient as a covariate in the model. The benefit of doing so is that unbiased variance components may be obtained because the variance component estimates may be inflated when inbreeding or heterozygosity are not fitted within the models [17,18,19]. According to Bolormaa et al. [4] and Xiang et al. [19], inbreeding or heterozygosity are fitted within the models to account for directional dominance (i.e., a higher percentage of positive than negative dominance effects), otherwise, the variance components for the dominance effects will be biased. Vitezica et al. [17] also reported that in the absence of inbreeding within the models, the estimates for dominance variance were inflated for litter size in a purebred pig line. In the present study, there was a slight decrease in the additive genetic variance when the models adjusted for inbreeding depression, suggesting that not fitting inbreeding in the model might generate overestimated additive genetic variances.

The heritability estimates obtained from the models including only the additive and dominance genetic effects (MAD) corroborate with the estimates reported in previous studies using the same public dataset [32]. The studies reported dominance variance aiming to propose different methods for obtaining additive and dominance genetic variance components [33,34,35]. The dominance variance ratio estimates for T3 reported by Da et al. [33] and Liu et al. [35] were greater than estimated in this study (0.07 vs. 0.02), and the dominance variance ratio estimated for T5 was smaller than that one reported by Nishio and Satoh [34] (0.045 vs. 0.063). Da et al. [36] proposed multifactorial methods with SNP and haplotypes to calculate the genomic relationship matrix and fitted the epistasis effects up to third-order in the genomic prediction models. When comparing the estimates from the most complete models, in general, Da et al. [36] reported similar heritability (0.02 vs. 0.01, 0.22 vs. 0.21, 0.13 vs. 0.12, 0.33 vs. 0.31, 0.36 vs. 0.34 for T1, T2, T3, T4, and T5, respectively) and lower additive-by-additive epistatic variance ratio estimates (0.04 vs. 0.13, 0.18 vs. 0.24, 0.28 vs. 0.29, 0.02 vs. 0.09, 0.05 vs. 0.13 for T1, T2, T3, T4, and T5, respectively). The differences between estimates reported by each study using the same dataset might be due to different methods and/or parametrization considered to obtain the genomic relationship matrices in each study.

The method to obtain the relationship matrices applied in this study, proposed by Vitezica et al. [6], is a flexible approach that can also be applied for populations that are not in Hardy-Weinberg equilibrium. However, it assumes linkage equilibrium between the genetic markers to ensure orthogonality between variance components. The non-orthogonality is clearly observed when there are significant changes in the variance estimates and when new variance components are introduced into the model [6]. In this study, the non-orthogonality between the variance components was noticeable when epistasis, especially additive-by-additive, was added to the model. This source of non-orthogonality might be due to linkage disequilibrium (LD) between the markers since this first dataset is from a nucleus purebred pig line, which is under intense selection pressure for various traits. Vitezica et al. [6] also reported high additive-by-additive variance and a sum of non-additive variances greater than the additive variance for a simulated high LD scenario due to strong selection. The additive, dominance, and epistasis effects are assumed independent of each other under the linkage equilibrium assumption [37, 38]. However, the presence of LD might affect the magnitude of dominance and epistasis variance due to the correlation between markers by introducing dependency between the genetic additive and non-additive values, which could lead to overestimation of the non-additive genetic variance components. As a result, the partition between additive and non-additive components is more difficult when LD between markers is present [39]. It is worth noting that many complex traits are predominantly additive [40], and the additive genetic variance often comprises most of the total genetic variance for most traits [41]. Accurate estimation of additive genetic variance components can be achieved without assuming linkage equilibrium. Therefore, incorporating non-additive genetic effects into genomic prediction may not offer significant advantages for selection purposes. However, to evaluate the genetic architecture of complex traits, there is still a need to develop methods for accurately estimating non-additive genetic variance components without assuming linkage equilibrium, especially considering that many livestock populations are subjected to intensive selective breeding (a major cause of LD).

As summarized by Wade et al. [42], for two bi-allelic loci, epistasis can be considered in different ways: (i) additive-by-additive, interactions between homozygotes at both loci; (ii) additive-by-dominance, interactions between homozygotes at locus A and heterozygotes at locus B and vice-versa; (iii) dominance-by-dominance, interactions between heterozygotes at both loci. The first dataset used in this study belongs to a purebred population, which could be one of the reasons why the variance components for additive-by-dominance and dominance-by-dominance epistasis are close to zero. Purebred populations tend to exhibit reduced heterozygosity, which can make it more challenging to accurately estimate additive-by-dominance and dominance-by-dominance epistatic effects [43]. The results indicate that although there may be those types of interactions between loci in the purebred population evaluated, they might be mostly captured by the additive genetic variance [41].

All models presented dispersion close to one, indicating that there is no inflation or deflation for the genomic estimated breeding values (GEBVs) [30]. The GEBVs also presented bias close to zero in all models, which is desirable, since biased estimates compromise the estimates of genetic trends and genetic gains in breeding programs [44]. However, the accuracy differed across traits and traits with higher (additive) heritability yielded greater accuracies, which agrees with the literature [45,46,47]. Different accuracies were estimated for each model, and there was a decrease in the accuracy when epistasis effects were incorporated in the models, which can be associated with the decrease in the additive genetic variance components for those models, since this variance component are used to compute the accuracy [44].

Although models including epistatic effects were significantly better (based on LRT parameter) for some traits (T2, T3, and T5), including epistasis had lower GEBV accuracies for traits with lower heritability (T1, T2, and T3). It is noticeable that models including epistatic effects had greater standard errors for the additive genetic variance and heritability estimates for T1, which is the trait with the lowest heritability. This suggests that the available data might not be sufficient to effectively distinguish between additive genetic and epistasis effects. Therefore, larger datasets may be required to obtain more accurate estimates of non-additive genetic effects [22], particularly for traits with low heritability.

Another factor influencing this inferior accuracy for models including epistasis may be due to the non-orthogonality between the components, especially between additive genetic effects and additive-by-additive epistasis. At the same time, when only additive genetic effects were fitted in the models, part of the additive-by-additive epistasis was captured by the additive genetic component, as previously reported in other studies [3, 48]. Even though models including additive-by-additive epistasis was significantly better (based on LRT parameter) and this effect accounted for a relevant proportion of phenotypic variance for some traits in the purebred population, including epistasis in the model may capture part of the additive genetic variance [3]. The ranking of selection candidates also changed for models as a consequence of the change in variance partition due to the inclusion of non-additive genetic effects in the models.

In general, models including epistasis presented the largest changes in ranking of the animals in comparison to the additive model. According to Piccoli et al. [49], the degree of reranking is directly associated with the accuracies of estimates and, at the same time, it is difficult to compare the reranking considering real datasets since the true breeding values are unknown. However, it is noticeable that epistasis plays an important role in genomic prediction and might change the selected animals when considered in a genetic evaluation program, which would impact genetic progress for the traits evaluated. There are studies investigating the impact of incorporating epistasis into genomic predictions aiming to develop methods that effectively integrate this effect into genomic prediction models (e.g., [50, 51]). Previous studies have shown that including all possible marker interactions in the prediction model does not improve the prediction ability [52, 53]. However, there is an improvement in prediction ability when only some interactions between markers are included in the model, being selected only those interactions with greater effects [52, 53].

Dataset 2: crossbred pig population

Since all heat stress-related traits had substantial additive genetic variability, including them in a selection program is feasible and can lead to genetic progress. Although, it should be noted that phenotypes measured in crossbred animals are expected to be more influenced by non-additive genetic effects than in purebred populations [21, 22]. Small non-additive genetic variances were observed for most of the traits evaluated in the crossbred dataset. Additionally, the estimates for dominance and epistasis variances had large standard errors which illustrates the difficulty of obtaining good estimates [17] and suggesting that the dataset was not large enough [19].

Besides HD and PS having considerable epistasis estimates, models including epistatic effects were significant better (based on LRT parameter) only for PS and no significant dominance variance was observed for any trait related to heat stress in the crossbred population. However, HD and PS are categorical traits and, threshold models [54,55,56] should be more suitable for genetically evaluating them. As highlighted by Alves et al. [3], there might be some limitations in estimating non-additive genetic effects for categorical traits.

Estimating non-additive genetic variances is also difficult, since they are often, or at least partially, confounded with other effects, such as common environmental or maternal effects [4]. In this study, it is noticeable that, for traits related to heat stress, there was a confoundment between permanent environmental and non-additive genetic effects. This was observed through a reduction in the variance of permanent environmental effects when non-additive genetic effects were included in the model. This observation, particularly regarding PS, is similar to the results reported by Vitezica et al. [17]. The authors also reported that including non-additive genetic effects in the model did not have a large impact on residual variances and the non-additive genetic effects were captured by the permanent environmental effects [17].

The findings from this study suggest that the gene actions for heat stress-related traits in this population are mainly additive, or at least most of the non-additive genetic effects are captured by the additive genetic component. Although the dataset used in this study came only from crossbred animals, the pure breeds used for crossbreeding are both considered maternal breeds (Large White and Landrace) and genetically related [57,58,59]. Thus, in crossings between lines or breeds with high genetic distance it would be expected higher heterosis and greater non-additive genetic variance estimates [21] in comparison to the present study.

Implications and future studies

Quantitative traits are controlled by many genes with additive and potentially non-additive gene actions, and their phenotypic expression are the result of many developmental and biochemical pathways comprised of loci networks that interact at the genetic and molecular levels, generating the epistasis effects [60]. Thus, the magnitude of the epistasis effects depends on how many pathways are involved in the expression of certain phenotypes and how these pathways are connected and interact between them. In the other hand, the magnitude of dominance effects depends on the deviations of genotypic values from breeding values for each locus. If the proportion of positive dominance effects is greater than the proportion of negative dominance effects, it is referred to as directional dominance [1, 9]. Heterosis and inbreeding depression are directly dependent on non-additive genetic effects [9, 17]. Although it is expected that non-additive genetic variation would be larger in crossbred populations [5, 9, 22, 61], most of the traits measured in the crossbred population did not present substantial non-additive genetic variance ratios. Additionally, fitting non-additive genetic effects in the models did not improve the prediction ability of breeding values for the purebred population. However, including non-additive genetic effects, especially additive-by-additive epistasis, changed the animal’s ranking and can therefore affect selection decisions.

Many studies have suggested that non-additive genetic effects are greater in traits related to fitness and adaptation than in morphological traits [62,63,64,65]. Considering the importance of selecting more climatically resilient animals, it is necessary to elucidate the genetic background of traits related to heat stress, especially for crossbred pigs. In our study, although models including non-additive genetic effects were not significant (based on LRT parameter) for most of the traits related to heat stress response, these effects might play an important role in gene action for those traits. Thus, it would be important to perform studies considering a larger dataset including records from purebred and crossbred populations and evaluate different crossings as well.

Conclusions

Including dominance effects in genomic prediction models did not improve the predictive ability of the models and most traits had none to low dominance variance. In the purebred population, low dominance variance ratios were observed, potentially due to reduced heterozygosity. Although models including additive-by-additive epistasis effects in genomic prediction exhibited significantly better fit for some traits, it led to lower accuracy when genomically predicting breeding values for traits with low heritability. For the evaluated populations and models, the inclusion of non-additive genetic effects did not improve genomic prediction accuracy. The non-orthogonality between variance components, especially between additive and additive-by-additive epistasis, suggests the presence of linkage disequilibrium challenges in the partition of additive and non-additive genetic effects. Therefore, there is still a need for developing methods able to accurately estimate non-additive variance components without assuming linkage equilibrium.

Traits related to heat stress in the crossbred population had substantial additive genetic variability, indicating their suitability for inclusion in selection programs. The difficulty in estimating non-additive genetic effects and the confounding with other effects, such as permanent environmental effects, is a challenge to obtain accurate estimates. Models incorporating non-additive genetic effects in genomic prediction for heat stress-related traits did not perform better than the additive model for most of the traits evaluated. Future studies should use larger datasets when they become available, including both purebred and crossbred animals, to better understand the contribution of non-additive genetic effects, particularly for traits related to adaptation, welfare, and resilience.

Methods

Dataset 1: purebred pig population

Purebred dataset description

We used a public dataset [32] including pedigree, genotypic, and phenotypic information from a single nucleus pig line from Pig Improvement Company (PIC, Hendersonville, TN, USA). This dataset consisted of 3,534 individuals with genotypes for 52,843 SNPs and phenotypes for five performance traits that represent a small number of phenotypes that are routinely collected from birth in the genetic nucleus and with heritability estimates ranging from 0.07 to 0.62. The phenotype was either pre-corrected for environmental factors and rescaled by correcting for the overall mean (traits 3, 4, and 5) or was a rescaled, weighted mean of corrected progeny phenotypes (traits 1 and 2). The data also includes a pedigree file with parents and grandparents of the genotyped animals [32]. The descriptive statistics of the traits are presented in Table 7.

Table 7 Descriptive statistics of a public dataset from a purebred pig population [32]

The quality control of genotype data consisted of removing SNPs with call rate lower than 0.90, minor allele frequency (MAF) less than 0.01, and extreme deviation from Hardy-Weinberg equilibrium (p-value < 10⁻⁵). After quality control, 40,828 SNPs and 3,534 individuals remained for further analyses. No animals were removed during the quality control due their genotypes were imputed [32], then all animals presented higher call rate (> 0.90).

Variance components estimation

Different models including or not inbreeding and non-additive genetic effect of dominance and epistasis were fitted to estimate variance components:

• MA: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{Za}+\varvec{\upepsilon }$;

• MAI: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\varvec{\upepsilon }$;

• MAE: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{Za}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\varvec{\upepsilon }$;

• MAIE: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\varvec{\upepsilon }$;

• MAD: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{Za}+\textbf{Zd}+\varvec{\upepsilon }$;

• MAID: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Zd}+\varvec{\upepsilon }$;

• MADE1: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{Za}+\textbf{Zd}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\varvec{\upepsilon }$;

• MAIDE1: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Zd}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\varvec{\upepsilon }$;

• MADE2: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{Za}+\textbf{Zd}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\textbf{Z}{\textbf{e}}_{\textbf{ad}}+\varvec{\upepsilon }$;

• MAIDE2: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Zd}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\textbf{Z}{\textbf{e}}_{\textbf{ad}}+\varvec{\upepsilon }$;

• MADE3: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{Za}+\textbf{Zd}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\textbf{Z}{\textbf{e}}_{\textbf{ad}}+\textbf{Z}{\textbf{e}}_{\textbf{dd}}+\varvec{\upepsilon }$;

• MAIDE3: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Zd}+\textbf{Z}{\textbf{e}}_{\textbf{aa}}+\textbf{Z}{\textbf{e}}_{\textbf{ad}}+\textbf{Z}{\textbf{e}}_{\textbf{dd}}+\varvec{\upepsilon }$;

where $\varvec{y}$ is the vector of observations; $\varvec{\upbeta }$ is the vector of fixed effects; $\textbf{f}$ is the vector of inbreeding coefficient based on pedigree; $\textbf{b}$ is the inbreeding depression; a is the vector of additive genetic effects with $\textbf{a} \sim N(0, {\textbf{G}}_{\textbf{A}}{\sigma }_{a}^{2})$, where ${\sigma }_{a}^{2}$ is the additive genetic variance and ${\textbf{G}}_{\textbf{A}}$ is the genomic additive relationship matrix; $\textbf{d}$ is the vector of dominance effects with $\textbf{d} \sim N(0, {\textbf{G}}_{\textbf{D}}{\sigma }_{d}^{2})$ where ${\sigma }_{d}^{2}$ is the dominance variance and ${\textbf{G}}_{\textbf{D}}$ the genomic dominance relationship matrix; ${\textbf{e}}_{\textbf{aa}}$ is the vector of additive-by-additive epistatic effects, ${\textbf{e}}_{\textbf{aa}} \sim N(0, {\textbf{G}}_{\textbf{AA}}{\sigma }_{aa}^{2})$, where ${\sigma }_{aa}^{2}$ is the additive-by-additive epistatic variance and ${\textbf{G}}_{\textbf{AA}}$ the genomic additive-by-additive epistatic relationship matrix; ${\textbf{e}}_{\textbf{ad}}$ is the vector of additive-by-dominance epistatic effects, ${\textbf{e}}_{\textbf{ad}} \sim N(0, {\textbf{G}}_{\textbf{AD}}{\sigma }_{ad}^{2})$, where ${\sigma }_{ad}^{2}$ is the additive-by-dominance epistatic variance and ${\textbf{G}}_{\textbf{AD}}$ the genomic additive-by-dominance epistatic relationship matrix; ${\textbf{e}}_{\textbf{dd}}$ is the vector of dominance-by-dominance epistatic effects, ${\textbf{e}}_{\textbf{dd}} \sim N(0, {\textbf{G}}_{\textbf{DD}}{\sigma }_{dd}^{2})$, where ${\sigma }_{dd}^{2}$ is the dominance-by-dominance epistatic variance and ${\textbf{G}}_{\textbf{DD}}$ the genomic dominance-by-dominance relationship matrix; $\varvec{\upepsilon }$ is the vector of residuals, $\varvec{\upepsilon } \sim N(0, \textbf{I}{\sigma }_{\epsilon }^{2})$, where ${\sigma }_{\epsilon }^{2}$ is the residual variance and $\textbf{I}$ an identity matrix; $\textbf{X}$ and $\textbf{Z}$ are incidence matrices of fixed and genetic effects, respectively.

The genomic additive relationship matrix (${\textbf{G}}_{\textbf{A}}$) was computed according to the method proposed by VanRaden [66]: ${\textbf{G}}_{\textbf{A}}=\frac{\textbf{MM}^{\prime}}{2\sum {\text{p}}_{\rm{i}}{\text{q}}_{\rm{i}}}$, where the M matrix contains elements equal to ($2-2{p}_{i}$), ($1-2{p}_{i}$), ($-2{p}_{i}$) for ${A}_{1}{A}_{1}$, ${A}_{1}{A}_{2,}$ and ${A}_{2}{A}_{2}$ genotypes, respectively, being ${p}_{i}$ the frequency of the allele ${A}_{1}$ in the SNP i and ${q}_{i}$ is equal to $1-{p}_{i}$. The genomic dominance relationship matrix (${\textbf{G}}_{\textbf{D}}$) was computed according to the method proposed by Vitezica et al. [5]: ${\textbf{G}}_{\textbf{D}}=\frac{\textbf{WW}'}{4\sum {p}_{i}^{2}{q}_{i}^{2}}$, where the W matrix contains elements equal to (-$2{q}_{i}^{2}$), ($2{p}_{i}{q}_{i}$), ($-2{p}_{i}^{2}$) for ${A}_{1}{A}_{1}$, ${A}_{1}{A}_{2}$, and ${A}_{2}{A}_{2}$ genotypes, respectively. The genomic epistatic relationship matrices were computed according to the methods proposed by Vitezica et al. [6]: ${\textbf{G}}_{\textbf{AA}}=\frac{{\text{G}}_{\text{A}}\odot\,{\text{G}}_{\text{A}}}{\text{tr}({\text{G}}_{\text{A}}\odot\,{\text{G}}_{\text{A}})/\text{n}}$, ${\textbf{G}}_{\textbf{AD}}=\frac{{\text{G}}_{\text{A}}\odot\,{\text{G}}_{\text{D}}}{\text{tr}({\text{G}}_{\text{A}}\odot\,{\text{G}}_{\text{D}})/\text{n}}$, and ${\textbf{G}}_{\textbf{DD}}=\frac{{\text{G}}_{\text{D}}\odot\,{\text{G}}_{\text{D}}}{\text{t}\text{r}({\text{G}}_{\text{D}}\odot\,{\text{G}}_{\text{D}})/\text{n}}$ being $\odot$ the Hadamard product, $tr$ the matrix trace, and $n$ the number of animals.

The analyses were performed using the ASREML software [67]. To understand how the phenotypic variance was split when different non-additive genetic effects were included in the models, variance components were estimated for each model. Estimates of additive heritability or narrow-sense heritability (${h}_{a}^{2}$), dominance variance ratio (${h}_{d}^{2}$), epistatic additive-by-additive variance ratio (${h}_{aa}^{2}$), epistatic additive-by-dominance variance ratio (${h}_{ad}^{2}$), and epistatic dominance-by-dominance variance ratio (${h}_{dd}^{2}$) were computed as the ratio of estimates of additive genetic variance ($\widehat{{\sigma }_{a}^{2}}$), dominance variance ($\widehat{{\sigma }_{d}^{2}}$), additive-by-additive epistatic variance ($\widehat{{\sigma }_{aa}^{2}}$), additive-by-dominance epistatic variance ($\widehat{{\sigma }_{ad}^{2}}$), and dominance-by-dominance epistatic variance ($\widehat{{\sigma }_{dd}^{2}}$) to the total phenotypic variance ($\widehat{{\sigma }_{p}^{2}}$).

Predictive ability

Bias, dispersion, and accuracy estimates were obtained using the Linear Regression (LR) method [44]. LR method is based on the comparison of breeding values estimated from partial datasets with a dataset containing all information (whole dataset). The whole dataset was composed of all animals and their observations, while ten partial datasets were generated by setting 10% of records missing to perform a 10-fold cross-validation. The prediction bias was computed as:

$${\mu }_{wp} = \overset{-}{{\widehat{\boldsymbol{u}}}_{p}} - \overset{-}{{\widehat{\boldsymbol{u}}}_{w}},$$

where ${\mu }_{wp}$ is the bias, $\overset{-}{{\widehat{\varvec{u}}}_{p}}$ is the GEBV mean predicted with partial information (GEBV_p) and $\overset{-}{{\widehat{\varvec{u}}}_{w}}$ is the GEBV mean predicted with whole information (GEBV_w). The GEBV dispersion is the slope of the regression (${b}_{w,p}$) of GEBV_w (${\widehat{\varvec{u}}}_{w}$) on GEBV_p (${\widehat{\varvec{u}}}_{p}$), and was obtained as:

$${b}_{w,p}=\frac{cov({\widehat{\varvec{u}}}_{p},{\widehat{\varvec{u}}}_{w})}{var\left({\widehat{\varvec{u}}}_{p}\right)}$$

The accuracy was obtained as:

$${acc}_{p}=\sqrt{\frac{cov({\widehat{\varvec{u}}}_{p},{\widehat{\varvec{u}}}_{w})}{\left(1-\overset{-}{F}\right)\widehat{{\sigma }_{a}^{2}}}}$$

where ${acc}_{p}$ is the accuracy, $\overset{-}{F}$ is the mean inbreeding coefficient based on pedigree from the animals with partial information, and $\widehat{{\sigma }_{a}^{2}}$ is the estimated additive genetic variance using the whole population.

Model comparison

The AIC [68] LRT parameters were used to compare the models. LRT was computed separately for the models that included or not the inbreeding coefficient as a covariate. This was done to specifically assess the impact of incorporating random non-additive genetic effects on the model, while isolating the influence of the inbreeding coefficient. The LRT statistic was calculated as:

LRT = − 2(LogL _{reduced model} − LogL _{complete model}),

in which LogL_{reduced model} and LogL_{complete model} are the likelihood logarithms for the reduced (i.e., from additive model; MA or MAI) and complete models (i.e., from the models including non-additive genetic effects), respectively. The level of significance used was 0.05.

We also evaluated the impact of the model used on the ranking of the animals and the proportion of commonly selected individuals based on the GEBVs from different models. The Spearman’s rank correlation was calculated to evaluate the re-ranking between the GEBVs from the MA model, which is the most used approach for GEBV calculation, and the other evaluated models. Assuming different selection intensities, the top 1%, 5%, 10%, and 20% of the animals were ranked according to their GEBVs based on the MA model and compared with the animals selected based on the models including non-additive effects.

Dataset 2: crossbred pig population

Crossbred dataset description

To evaluate the contribution of non-additive genetic effects on traits related to heat stress in lactating sows, we used a dataset previously described by Johnson et al. [68]. In summary, a total of 1,645 multiparous lactating sows (Large White × Landrace) were housed in either a naturally ventilated or mechanically ventilated farrowing barn. Thermoregulatory and anatomical traits were measured in 1,381 sows. The traits included in this study were skin surface temperature from ear (T_ES), shoulder (T_SS), rump (T_RS), and tail (T_TS); vaginal temperature considering the whole data measured every 10 min (T_Vall), and the average of the six records per hour corresponding to 08:00, 12:00, 16:00, and 20:00 h during four days (T_V4days); respiration rate (RR), panting score (PS; score scale from 0 to 3), and hair density (HD, score scale from 0 to 2) (the data collection is described in details by Johnson et al. [69]). All animals were genotyped using the PorcineSNP50K (50,703 SNPs) Bead Chip (Illumina, San Diego, CA, USA). Quality control of genotype data consisted of removing SNPs with a call rate lower than 0.90 and minor allele frequency (MAF) less than 0.01. It was also removed animals with call rate lower than 0.90. After quality control, 49,536 SNPs for 1,625 animals remained for further analyses. The descriptive statistics for the continuous traits are shown in Table 8.

Table 8 Descriptive statistics for continuous heat stress indicators in Landrace x Large White crossbred pig population

Estimation of variance components for traits related to heat stress

As performed for Dataset 1, we also fitted different models including or not non-additive genetic effect of dominance and epistasis to estimate variance components. However, we did not include additive-by-dominance and dominance-by-dominance epistasis effects to avoid overparametrized models, since this second dataset is small.

For the trait HD, the models MAI, MAID, MAIDE1, and MAIE previously mentioned were used. For skin temperatures, vaginal temperatures, PS, and RR repeatability models were fitted:

• MAIpe: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Zpe}+\varvec{\upepsilon }$;

• MAIDpe: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Zd}+\textbf{Zpe}+\varvec{\upepsilon }$;

• MAIDEpe: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Zd}+\textbf{Z}{\varvec{e}}_{\varvec{a}\varvec{a}}+\textbf{Zpe}+\varvec{\upepsilon }$;

• MAIEpe: $\textbf{y}=\textbf{X}\varvec{\upbeta }+\textbf{fb}+\textbf{Za}+\textbf{Z}{\varvec{e}}_{\varvec{a}\varvec{a}}+\textbf{Zpe}+\varvec{\upepsilon }$;

where $\textbf{f}$ is the vector of genomic inbreeding coefficient; $\textbf{b}$ is the inbreeding depression; $\varvec{p}\varvec{e} \sim N(0, \varvec{I}{\sigma }_{pe}^{2})$ where ${\sigma }_{pe}^{2}$ is the permanent environmental variance and $\varvec{I}$ an identity matrix; all other effects and matrices in the models have been previously described. The fixed effects for each trait were included as previously defined by Freitas et al. [70]. For T_Vall, the fixed effects were the concatenation of week and day of measurement; parity; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. For T_V4days, the fixed effects were concatenation of week, day, and time of measurement; parity; days in lactation; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. The fixed effects fitted for skin surface temperature (T_ES, T_SS, T_RS, and T_TS) and RR were trait recorder; concatenation of week, day, and time of measurement; parity; days in lactation; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. For PS, the fixed effects were trait recorder; concatenation of week and day of measurement; parity; days in lactation; concatenation of barn type and room; and in-barn environmental temperature as a linear covariate. The fixed effects fitted for HD were trait recorder and parity.

The variance components for dataset 2 were estimated following the same methods described for dataset 1. The analyses were also performed using ASREML software [67]. The estimates of additive heritability (${h}_{a}^{2}$), dominance variance ratio (${h}_{d}^{2}$), and epistatic additive-by-additive variance ratio (${h}_{aa}^{2}$) were computed, as described for dataset 1. The ratio of the total genetic variance ($\widehat{{\sigma }_{g}^{2}}=\widehat{{\sigma }_{a}^{2}}+\widehat{{\sigma }_{d}^{2}}+\widehat{{\sigma }_{aa}^{2}}$) explained by dominance ($\widehat{{\sigma }_{d}^{2}}/\widehat{{\sigma }_{g}^{2}}$) and epistasis additive-by-additive ($\widehat{{\sigma }_{aa}^{2}}/\widehat{{\sigma }_{g}^{2}}$) were also estimated.

Model comparison

The AIC [68] and LRT parameters were also used to compare the models’ adjustment. LRT was computed comparing the models including non-additive genetic effects with the additive animal model (MAIpe). The statistic was calculated as previously described and the level of significance used was 0.05.