Datasets from five nucleus herds located in North Carolina and Texas (USA) were provided by Smithfield Premium Genetics (SPG; Rose Hill, NC, United States) and consisted of records from purebred Landrace pigs. Measurements were collected between the years 2014–2019. The five traits included in this study were categorized as: 1) growth: birth weight (BW) and weaning weight (WW); 2) body composition: off-test backfat (BF); and 3) reproduction traits: total number of piglets born (TNB) and number of piglets born alive (NBA).

Records outside of an interval of 3.5 standard deviations (SD) from the mean were removed from further analyses. Observations from animals above the sixth parity or sows with less than 3 or more than 20 piglets per litter were also removed due to lower incidence. Animals with phenotypic observations but not present in the pedigree file and animals whose pedigrees were not complete for at least four generations were also removed from the analyses. This quality control was done to keep only relevant information for the variance component estimation analyses (Nilforooshan and Saavedra-Jiménez, 2020), especially for estimating epigenetic effects that require deeper knowledge of the genealogy and more complex partition of the phenotypic variance (Tal et al., 2010). The “optiSel” package in R (Wellmann, 2019) was used to check for erroneous pedigree entries and pedigree completeness. After the phenotypic data editing, the contemporary groups (CG) were defined by the concatenation of birth year, month, and farm for the growth and body composition traits and farrow year, month, and farm for the reproductive traits. The CG were defined based on the trait groups to better capture animals that would move together and share similar management practices (i.e., contemporaries). The number of animals and the descriptive statistics for all traits, after quality control, as well as the number of CG are included in Table 1. A minimum number of five animals was set for each level of the CG effect as well as for all other fixed effects included in the models (described later).

After pedigree quality control, the pedigree consisted of 17,794 animals, spanning a maximum of 12 generations, with a mean (SD) pedigree completeness index of 0.84 (0.33) for growth and body composition traits, while for reproductive traits, the pedigree comprised 8,390 animals, spanning 11 generations, with a pedigree completeness index of 0.83 (0.29). The “optiSel” package in R (Wellmann, 2019) was also used to calculate the pedigree statistics shown above. The RENUMF90 package from the BLUPF90+ family of programs (Misztal et al., 2018) was used to renumber the pedigree and keep up to four generations back from the animals with phenotypic and/or genomic information to reduce the pedigree and, consequently, time to reach convergence and memory usage during the analyses (Nilforooshan and Saavedra-Jiménez, 2020). Previous analyses tracing up to 20 generations back in the pedigree using RENUMF90 (Misztal et al., 2018) did not affect the estimates (results not shown).

Genotypes were available for 12,759 individuals genotyped with a Porcine SNP50K Bead Chip (Illumina, San Diego, CA, United States). The quality control (QC) was performed using the PREGSF90 software (Misztal et al., 2018) to remove samples and SNPs with call rate lower than 0.90 and SNPs with minor allele frequency lower than 0.05, SNPs with a difference between expected and observed heterozygous greater than 0.15, and SNPs located in non-autosomal chromosomes or with unknown genomic positions. After QC, 12,759 individuals and 34,524 SNPs remained for further analyses.

The transgenerational epigenetic component of covariance between relatives was derived based on the method proposed by Tal et al. (2010). The main assumption made by the authors was that the covariances between close relatives change when epigenetic marks are being transmitted across generations. Later, Varona et al. (2015) used this theory to derive the Λ and Λ ⁻¹ . The Λ matrix is calculated using the auto-recursive parameter λ , in which λ = 0.5 1 − ν , where ν is the reset coefficient, defined as the probability that the epigenetic state of the genome changes from one generation to the next. In brief, Tal et al. (2010) proposed a reset coefficient ( ν ) of the epigenetic marks and its complement ( 1 − λ ) as the transmissibility of the epigenetic marks, which weights the covariance between relatives in the presence of factors that can cause epigenetic changes. A Python function (see Additional file 1) was implemented to calculate the lower triangular and diagonal portion of Λ ⁻¹ directly in the long format to estimate the transgenerational epigenetic component of variance ( σ ξ 2 ) and heritability ( h ξ 2 ) using the BLUPF90 family programs (Misztal et al., 2018). The Phyton function developed in this study was based on the algorithm proposed by Varona et al. (2015), which assumed the same epigenetic variance across generations and between sexes. These assumptions allow the calculation of the epigenetic heritability (or epigenetic ratio) as h ξ 2 = σ ξ 2 ^ σ y 2 ^ , where σ y 2 ^ is the phenotypic variance, i.e., is the sum of all variance components from the model used to evaluate each trait.

Model building was conducted before fitting the epigenetic models. The systematic effects (presented later) were variables statistically significant at 5% probability in a linear model from a set of environmental (non-genetic) variables. Once the systematic effects were defined for each trait, the set of random effects to be used in the genetic models were chosen under the Restricted Maximum Likelihood (REML) (Patterson and Thompson, 1971) approach using the Average Information algorithm (AI-REML) (Gilmour et al., 1995) and the Akaike Information Criterion (AIC) (Akaike, 1974). The rank for the genetic models using different sets of random effects is presented in Supplementary Material 2: Supplementary Table S1, and the best genetic model (lowest AIC with meaningful estimates of the genetic parameters, i.e., within or close to previous estimates for the parameter and trait) for the evaluated traits were: y = Xb + Zu + Zm + Sq + e (1) y = Xb + Zu + Zm + e (2) y = Xb + Zu + Wpe + e (3) where Equations 1–3 are the genetic models for the growth (BW and WW), back fat thickness (BF), and reproductive (TNB and NBA) traits, respectively; y is the vector of the phenotypic records for each trait; b is the vector of systematic effects (for BW: gender, birth parity, and CG; WW: gender, birth parity, CG, and weaning age and birth weight as covariates; BF: gender, birth parity, CG, and off-test age as covariate; TNB and NBA: farrowing age as linear and quadratic covariates, CG, and parity); u is the vector of random direct additive genetic effects; m is the vector of random maternal additive genetic effects; q is the vector of random common litter environment effects; p e is the vector of random permanent environmental effects; and e is the vector of random residuals. X , Z , S , and W are the incidence matrices for b , u , m , q , and p e , respectively. BW was used as a covariate in the WW to account for the starting point of the growth trajectory in this trait. The AIREMLF90 software (Misztal et al., 2018) was used to estimate variance components under the following assumptions for each genetic model: u m q e ∼ N 0 0 0 0 , A σ u 2 A σ u , m 0 0 A σ u , m A σ m 2 0 0 0 0 I σ q 2 0 0 0 0 I σ e 2 (4) u m q e ∼ N 0 0 0 0 , A σ u 2 0 0 0 0 A σ m 2 0 0 0 0 I σ q 2 0 0 0 0 I σ e 2 (5) u m e ∼ N 0 0 0 , A σ u 2 A σ u , m 0 A σ u , m I σ m 2 0 0 0 I σ e 2 (6) u p e e ∼ N 0 0 0 , A σ u 2 0 0 0 I σ pe 2 0 0 0 I σ e 2 (7) where the assumption (Equation 4) was used for BW, Equation 5 for WW, Equation 6 for BF, and Equation 7 for TNB and NBA. The pedigree-based relationship matrix ( A ) was used to model the covariances in u and m in the genetic analysis, and the identity ( I ) matrix was used for the q , p e , and e . The σ u 2 , σ m 2 , σ q 2 , σ pe 2 , and σ e 2 are the variance components for u , m , q , pe , and e , respectively, and σ u , m is the covariance between u and m. The genetic models, i.e., models using the A as covariance structure for the additive genetic effects, will be called BLUP models from now on.

The A matrix was replaced by the H matrix in the BLUP models when the genomic information was included in the estimation of genomic-based variance components, i.e., single-step genomic BLUP models (ssGBLUP) (Legarra et al., 2009; Christensen and Lund, 2010). This H is the matrix that combines the pedigree and genomic information, and its inverse ( H − 1 ) was computed directly as (Aguilar et al., 2010): H − 1 = A − 1 + 0 0 0 τ α G + β A 22 − 1 − ω A 22 − 1 (8) where A − 1 , A 22 , A 22 − 1 , and G are the inverse of the pedigree relationship matrix, part of the A related to the genotyped animals, its inverse, and the genomic relationship matrix, respectively (Equation 8). The G was computed as in method 1 proposed by VanRaden (2008). The scaling parameters ( τ and ω ) were equal to 1.00, and the blending parameters α and β were equal to 0.95 and 0.05, respectively.

The epigenetic model for each trait was obtained by expanding the BLUP models for each trait including the transgenerational epigenetic effect, which will be called Epi-BLUP. In summary, the term Z ξ was included in the models (Equations 1–3), which became: y = Xb + Zu + Zm + Sq + Z ξ + e (9) y = Xb + Zu + Zm + Z ξ + e (10) y = Xb + Zu + Wpe + Z ξ + e (11) where Equations 9–11 are the epigenetic models for growth (BW and WW), backfat thickness (BF), and reproductive (TNB and NBA) traits, respectively; ξ is the vector of epigenetic effects, with ξ   ∼ N 0 , Λ σ ξ 2 assumed to be uncorrelated to all the other random effects, as proposed by Tal et al. (2010) and Varona et al. (2015); all the other effects were previously described. To estimate the transgenerational epigenetic heritability, the best epigenetic model, i.e., the model that presented the best λ , was defined first because the reset and transmissibility of the epigenetic markers are population-specific parameters. In this sense, a grid search ranging from 0.05 to 0.45 by 0.05 (totaling nine epigenetic models for each trait) was evaluated for λ in the creation of Λ ⁻¹ that was fitted directly in the epigenetic models. The H matrix (Legarra et al., 2009) was also used to model the additive genetic relationships in the epigenetic models to study the impact of genomic information in these models, i.e., epigenetic models including genomic information (Epi-ssGBLUP) and not epigenomic models (e.g., using whole-genome DNA methylation).

Even though both REML-based and Bayesian inference-based methods should converge to the same population parameters, there are differences in the statistical properties of these methods, which are important to consider, especially if the models are complex (Gianola and Fernando, 1986). In this context, we have estimated variance components based on BLUP, ssGBLUP, Epi-BLUP, and Epi-ssGBLUP models using Bayesian inference. Bayesian inference was used in this step because the epigenetic models have a complex covariance structure and this method is more robust to model complexity (Gianola and Fernando, 1986).

The best Epi-BLUP or Epi-ssGBLUP model for each trait, i.e., with the best λ , under Bayesian inference was chosen using the Deviance Information Criterion (DIC; lower values indicate better model fit) (Spiegelhalter et al., 2002). Maximum likelihood estimates from the AI-REML runs were used as starting values in the Bayesian inference to improve the sampling process (Harms and Roebroeck, 2018) as poor starting values tend to slow model convergence (Raftery and Lewis, 1992). The THRGIBBS1F90 software (Misztal et al., 2018) was used to estimate variance components considering all traits as linear traits. Initial values for size, burnin, and interval in the MCMC chain were defined as 100,000, 10,000, and 10, respectively, for all traits and models, and these values were increased if convergence was not achieved based on the author’s experience and final values are presented in the results section. Convergence was verified based on the Geweke diagnostic (Geweke, 1992) with 5% probability, a stationarity test (Heidelberger and Welch, 1983) criterion, and visual inspection. The convergence tests and visual inspection of the chains were done using the R software (R Language, 2012) with the “boa” package (Smith, 2007).

After estimating the variance components, the solutions of the mixed model equations for the best fit models were evaluated. The purposes of these analyses were mainly to: 1) investigate the association between the solutions of the additive, maternal, and permanent environment effects in a genetic or genomic model including the transgenerational epigenetic effects; 2) assess the changes in the solutions by including the transgenerational epigenetic effects; and 3) evaluate the prediction accuracy, bias, and dispersion of the breeding values in young non-phenotyped individuals when including the transgenerational epigenetic effects in the models.

The association between the additive, maternal, and permanent environment effects when including the transgenerational epigenetic effects in the models was evaluated using a Pearson correlation between the solutions considering all animals used for variance component estimation (please see the sections Traits and dataset edits and Pedigree and Genotypes). Changes in the solutions when including the transgenerational epigenetic effects were evaluated using a paired t-test with a significance level of 0.05% for the difference between the solutions, also considering all animals used for variance component estimation.

The prediction accuracy, bias, and dispersion were investigated using the Linear Regression method (LR) (Legarra and Reverter, 2018). In brief, a set of animals born in 2019 (N = 571) chosen among the ones used for the estimation of variance components had their phenotypic records masked for BW, WW, and BF, and were considered as the focal animals (young non-phenotyped selection candidates). In the case of TNB and NBA, the focal animals for the LR method were born in 2017 (N = 933), because this was the last year with phenotyped animals for the reproductive traits. After defining the focal animals, the whole and partial datasets for applying the LR method were created. The whole dataset included all animals and phenotypes used for variance component estimation in the prediction, while in the partial dataset the phenotypes for the focal animals were removed, so that their estimated breeding values (EBV) were obtained based on the relationships with the remaining animals with phenotypic records. In the end, the prediction accuracy, bias, and dispersion of the EBV for the focal animals were calculated as: Accuracy = cov EB V W , EB V I 1 − F ¯ σ ^ u 2 (12) Bias = ave u EB V I ^ − ave u EB V W ^ (13) Dispersion = cov EB V W , EB V I var EB V I − 1 (14) where cov EB V W , EB V I is the covariance between the EBV in the whole ( EB V W ) and incomplete, i.e., partial, ( EB V I ) datasets, F ¯ is the average inbreeding (pedigree or genomic), ave represent the arithmetic average function, u EB V W ^ and u EB V I ^ are the predicted EBV in the whole and incomplete datasets, respectively, and v ar EB V I is the variance of the EBV in the incomplete datasets (Equations 12–14). The other equation components were previously described.

The Epi-BLUP model presented a better fit (lower DIC) compared to BLUP for BW, WW, and BF (Table 2). λ = 0.05 provided the lowest DIC for BW and WW in the Epi-BLUP models (10,387.47 and 44,350.08, respectively), and λ = 0.10 was the best for BF (−30,458.31), which were lower than the DIC of the BLUP models for these traits (21,258.82, 52,803.60, and −17,142.00, respectively). Epi-BLUP models with λ = 0.05 provided slightly lower DIC for TNB and NBA (52,651.09 and 51,609.69, respectively) compared to the BLUP model (52,651.31 and 51,610.34, respectively), suggesting that the Epi-BLUP models would be a better choice for these traits as well.

When including genomic information, DIC values were lower in the best Epi-ssGBLUP models for BW, WW, BF, TNB, and NBA (−1,208,421.73, −924,295.31, −6,259,876.37, −112,151.59, −105,368.96, respectively) compared to the ssGBLUP models for these traits (−453,655.16, −439,644.93, −377,412.68, −11,668.76, −101,485.82, respectively). Similar ranks were observed between the Epi-ssGBLUP and Epi-BLUP models under Bayesian inference, with the best λ value of 0.05 for BW, WW, TNB, and NBA, except for BF, with the best λ value of 0.15.

The results of the genetic parameter estimation for BW are presented in Table 3. Chain size was much higher in the epigenetic models to obtain stable samples in the Markov-Monte Carlo chain (MCMC). The size of the chain for the Epi-BLUP model (biggest chain) was 1,750,000 samples with a burn-in of 1,000,000 and thin of 100, while for the BLUP model (smallest chain), these parameters were 220,000, 10,000, and 10, respectively. Small to no differences were observed among the variance components and genetic parameters between BLUP, Epi-BLUP, and ssGBLUP for BW, except for σ e 2 ^ , which was much smaller in the Epi-BLUP model. In this sense, σ P 2 ^ was also similar across the BLUP, Epi-BLUP, and ssGBLUP models, indicating that there was a partition of the variances in the Epi-BLUP model which resulted in the σ ξ 2 ^ being extracted directly from the residuals. In the case of the Epi-ssGBLUP, there was an inflation of the σ u 2 ^ and σ m 2 ^ compared to the other models, which were more than two and three times larger, respectively, in the Epi-ssGBLUP. Despite the reduction in the proportion of the phenotypic variation explained by the common litter effect ( q 2 ) and σ ξ 2 ^ in the Epi-ssGBLUP model, the σ P 2 ^ in this model was the largest compared to all other models for BW (0.766 vs. ∼0.560). The additive direct and maternal heritability estimates followed a similar pattern of the variance components across models, but the q 2 and σ ξ 2 , and consequently h ξ 2 , were higher in the Epi-BLUP model (0.296, 0.186, and 0.330, respectively) compared to the Epi-ssGBLUP model (0.204, 0.160, and 0.209, respectively). A similar pattern of the results was observed in the WW and BF (Tables 4, 5, respectively).

The results of the MCMC parameters, variance components, and genetic parameters for TNB and NBA are presented in Tables 6, 7, respectively. Similar patterns were observed for the MCMC parameters, σ u 2 ^ , σ P 2 ^ , h u 2 , and h ξ 2 in the TNB and NBA compared to the growth and body composition traits. However, the source of σ ξ 2 ^ in the Epi-BLUP and Epi-ssGLUP models was not due to σ e 2 ^ but due to the σ pe 2 ^ in the reproductive traits. There were no substantial differences in the σ e 2 ^ across all models for TNB and NBA. On the other hand, σ pe 2 ^ reduced in both Epi-BLUP and Epi-ssGLUP models by the same amount of variance captured by σ ξ 2 ^ in TNB and NBA. Also, similar to what was observed for growth and body composition traits, the σ P 2 ^ was similar to the BLUP, Epi-BLUP, and ssGBLUP estimates, suggesting that there was only repartitioning of variances in the Epi-BLUP compared to BLUP and ssGBLUP. The inflation in the σ u 2 ^ , σ P 2 ^ , and h u 2 estimates from the Epi-ssGBLUP for the reproductive traits was also noticeable, with 3-fold increases for σ u 2 ^ and h u 2 , and 14% for σ P 2 ^ in comparison to the ssGBLUP estimates. The posterior means of the h ξ 2 for TNB and NBA in the Epi-BLUP (0.044 and 0.042, respectively) and Epi-ssGBLUP (0.033 and 0.032, respectively) were much lower compared to those for BW, WW, and BF (0.330, 0.317, 0.336, respectively, for the Epi-BLUP and 0.436, 0.229, and 0.209 for the Epi-ssGBLUP, respectively).

Despite the assumptions of independence among the random effects in all models, we observed that the posterior samples are correlated in most cases (Figure 1). The magnitude of the correlations among the variance components tended to be smaller in the Epi-BLUP and Epi-ssGBLUP in comparison to the estimates from the BLUP and ssGBLUP models for all traits, respectively, except for σ ξ 2 , which was not fitted in these later models. Similar correlations were observed between the variance components in the BLUP and ssGBLUP models. The σ ξ 2 was highly negatively correlated to σ e 2 in the BW and WW for both Epi-BLUP and Epi-ssGBLUP models ( r ≤ −0.99), while only a low negative correlation was observed between σ ξ 2 and σ u 2 for these traits (−0.23 ≤ r ≤ −0.07). For BF, σ ξ 2 was also highly negatively correlated to σ e 2 with Epi-BLUP ( r = −0.98), while a high negative correlation was observed for Epi-ssGBLUP ( r = −0.80), and a negative moderate correlation was observed between σ ξ 2 and σ u 2 in the Epi-BLUP and Epi-ssGBLUP models (−0.34 and −0.55, respectively). A low negative correlation was observed between σ ξ 2 and σ e 2 for TNB and NBA in both Epi-BLUP and Epi-ssGBLUP models (−0.15 ≤ r ≤ −0.10), while a moderate to high negative correlation was observed between σ ξ 2 and σ pe 2 (−0.87 ≤ r ≤ −0.70), and a low to moderate negative correlation was observed for σ ξ 2 and σ u 2 (−0.32 ≤ r ≤ −0.18).

The correlations between solutions for u , m , and ξ for BW, WW, and BF, and u , pe , and ξ for TNB and NBA in the BLUP and Epi-BLUP models are presented in Figure 2. Only the correlations among solutions in the genetic and epigenetic models with the best λ , not including genomic information, are presented here. High correlations were observed between the direct additive genetic solutions in BLUP and Epi-BLUP models for all traits, in which the lowest correlation of 0.85 was observed for BF and all others were above 0.96. The correlation for the maternal genetic solutions between BLUP and Epi-BLUP models was low for BF (−0.22) and high for BW and WW (0.99 and 0.97, respectively). In the case of TNB and NBA, the permanent environment solutions were highly correlated (0.98) between BLUP and Epi-BLUP models.

A paired t-test at 5% significance level was applied to evaluate the mean difference among the solutions from the BLUP and Epi-BLUP models for the effects mentioned above and are presented in Figure 3. The mean difference between the direct additive and maternal genetic solutions between the BLUP and Epi-BLUP models for all traits were statistically significant. No statistical significance was observed for the mean difference between the permanent environmental solutions comparing the results from the BLUP and Epi-BLUP models.

Only the genetic and best epigenetic models were used for comparing the prediction results using the LR method. Prediction accuracies ranged between 0.199 and 0.443 for the Epi-BLUP with λ = 0.05 in the BW and Epi-BLUP with λ = 0.10 in the BF, respectively (Figure 4). Prediction bias ranged between −0.080 and 0.034 for the BLUP model on WW and NBA, respectively. The dispersion ranged between −0.134 and 0.131 for the WW BLUP and TNB Epi-BLUP with λ = 0.05, respectively. No substantial differences were observed between the BLUP and Epi-BLUP models for the traits evaluated based on all predictivity measures (accuracy, bias, and dispersion), indicating that the inclusion of transgenerational epigenetic effects did not affect the prediction of the breeding values in the focal animals.

Epigenetic models are quite complex. Even though there is the assumption that there is independence among genetic, epigenetic, and residuals (Tal et al., 2010; Varona et al., 2015), non-negligible linear association between those components exists, as shown in Figure 1. Varona et al. (2015) also reported a posterior association between the transgenerational epigenetic and the additive genetic and residual variance components. Even though there were previous studies that used REML-based methods to estimate transgenerational epigenetic variances (Paiva et al., 2018a; Paiva et al., 2018b), Bayesian inference was preferred to conduct the analysis in this study. Among the main properties of Bayesian inference, it can be highlighted the fact that the variance components are never out of parameter space (Box and Tiao, 1973), and inferences are exact for any sample size (Gianola and Fernando, 1986). We also tried to fit AI-REML-based epigenetic models, but the results were not stable. The number or iterations to converge was very large for some models/traits using AI-REML and many epigenetic models failed to converge (Supplementary Material S3), which is expected as REML estimation is sensitive to complex models (Gianola and Fernando, 1986; Misztal, 2008). The advantages of Bayesian- over REML-based methods make the former more robust for complicated models than the latter. In this sense, only Bayesian inference results were kept in the main text.

All epigenetic models, including or not genomic information, using Bayesian inference converged (Table 2) given the criteria used and led to similar best λ values across traits (0.05 for all traits, except BF, which was 0.10). Epigenetic reprogramming can be produced in laboratory in many ways (Cantone and Fisher, 2013; Basu and Tiwari, 2021), but naturally happens twice during mammalian life (gametogenesis and preimplantation development) (Lacal and Ventura, 2018). In the context of this study, we expect that epigenetic reprogramming is less likely to be influenced by traits such as litter size or backfat depth in sows; rather, it is more likely to influence these traits. Among the factors that affect epigenetic reprograming, there are reports of sex (Lacal and Ventura, 2018), age (Caulton et al., 2022), and environmental factors/stress (Tal et al., 2010).

Even though the Epi-ssGBLUP provided lower DIC, the Epi-BLUP models were chosen as the best because of the clear inflation in the additive genetic, hence phenotypic, variances in the epigenetic models including genomic information. Large changes in variance components previously estimated by introducing new random effects indicate non-orthogonality of the effects in the model (Vitezica et al., 2017), which is related to not accounting for significant interaction between random effects, i.e., their covariance matrices (Zhou et al., 2020). This was not observed in the pedigree-based epigenetic models for all traits (Tables 3–7), i.e., no substantial changes were observed in the variance components when adding the transgenerational epigenetic effect to the genetic model (moving from the BLUP to the Epi-BLUP model). Orthogonality is a very important and useful property because it directly affects the partition of the genetic effects (Vitezica et al., 2017), disqualifying epigenetic models including genomic information (Epi-ssGBLUP) for variance component estimation. Given the current derivation of the epigenetic relationship matrix (Tal et al., 2010; Varona et al., 2015) and the analysis framework employed, i.e., independence of genetic and epigenetic effects, the Epi-BLUP models are more reliable than Epi-ssGBLUP. As the variance components and genetic parameters using a genetic or genomic model (BLUP and ssGBLUP, respectively) were not substantially different, the BLUP model was used as the benchmark for comparisons with Epi-BLUP models because they use the same additive genetic relationship matrix. Similarities between the variance components and genetic parameters estimated using pedigree and genomic models were also reported in the literature by Aldridge et al. (2020), which recommended the use of the pedigree relationship matrix for variance component estimation due to easier implementation.

Modeling the covariance between epigenetic and other random effects in the model could be beneficial, especially for the Epi-ssGBLUP, where the non-orthogonality was evident. To model the interaction between additive genetic effects based on genomics and transcriptomic data, Perez et al. (2022) used Hadamard products between the covariance matrices of these effects and, as an alternative, developed a model where these effects were independent. The independence between the correlated genomic and transcriptomic information was achieved by decorrelating these effects using a “smoother matrix,” which relied on the assumption that the transcriptomic was conditioned to the genomic information, which was called the GTCBLUP model (Perez et al., 2022). Using an interaction effect in the model resulted in non-stable estimates of the variance components, while the GTCBLUP model resulted in a stable partition of the variance components (Perez et al., 2022). Since epigenetic variance originates from epigenetic marks themselves, applying a model similar to that of Perez et al. (2022) would require incorporating epigenomic information (e.g., DNA methylation) into the animal model and conducting the appropriate statistical analysis. Unfortunately, no epigenomic information was present for the animals included in the study. Zhou et al. (2020) proposed a model that considers correlations between random effects in genetic models and showed the importance of modeling the covariance between these effects when non-negligible associations exist (not always linear associations). The core of the work of Zhou et al. (2020) was in providing a framework for multiomics data analyses, but it can be applied in the future for fitting epigenetic models including genomic information. This is because, in the end, the problem of fitting these models is the same: similarity or dissimilarity between the covariance matrices associated with the random effects in the model. However, more theoretical work to understand the association between transgenerational epigenetic and genomic information is required to make assumptions about how to model these effects and their interactions, which were beyond the scope of the present work.

Epigenetics affects the phenotypic variation through changes in gene expression (Heard and Martienssen, 2014), which means that even though the epigenome is not expected to change the DNA sequence, it is modulating or interacting with it during phenotypic expression. In the same way, as epigenetic changes can be stress- or environmentally induced (Ibeagha-Awemu and Yu, 2013), non-negligible associations between epigenetic and residual effects could be produced. Just assuming independence between epigenetic and other effects affecting the phenotypic variation resulted in simpler mathematical work in deriving the covariance between relatives from a quantitative genetics perspective. Even though not optimal from a biological point of view, the assumptions originally made can be considered robust enough in epigenetic models not including genomic information as these models met desired expectations, e.g., orthogonality of the effects. The problems related to non-orthogonality between the additive genetic and epigenetic effects in the genomic analysis (between H and Λ instead of between A and Λ ) may be related to the fact that the epigenetic relationships were derived based on the pedigree-based average relationships. Genomic information provides realized genomic relationships at the QTL level, which changes the additive genetic relationship matrix by adjusting the average relationships to observed values and propagating this information to non-genotyped related to genotyped individuals (Legarra et al., 2009). Changes in the average relationships in closely related individuals plus adding relationship coefficients to animals without pedigree links in the H matrix may have broken important assumptions to hold non-orthogonality when genomic information was added to the epigenic models.

Assuming 0.05 to 0.10 and the best λ values (range of best λ across traits – Table 2), 10%–20% of the epigenetic markers are being inherited in the studied population, i.e., 80%–90% of the epigenetic marks are being reset. No previous reports for the λ , reset, or percentage of epigenetic marks inherited in pigs were found in the literature. The expectation is that the percentage of epigenetic marks being inherited should not be high because it results from an error in resetting the epigenetic marks from one generation to the next (Tal et al., 2010; Basu and Tiwari, 2021), which agrees with our results. Nevertheless, new studies including epigenetic marks (DNA methylations, chromatin remodeling, miRNA, among others) at the individual level could provide insights into the estimation of the reset and inheritance of epigenetic marks.

The variance components and genetic parameters from the BLUP and ssGBLUP models were consistent with the values observed in the literature for all traits evaluated. Alves et al. (2018) reported heritability estimates of 0.05 and 0.07 for the direct additive genetic effects in BW of Canadian Landrace pigs with two datasets. The maternal genetic heritability for the BW presented by Alves et al. (2018), 0.27, was higher than the one observed in this study (∼0.13 for the BLUP, ssGBLUP, and Epi-BLUP models – Table 3), while the proportion of the phenotypic variation explained by the common litter effect was lower (0.17, compared to ∼0.30 from the BLUP, ssGBLUP, and Epi-BLUP models). Even though the base genetic models were the same (fixed effects plus direct and maternal genetic additive and common litter effects), the population and data structure were different (e.g., different numbers of individuals with phenotypes and the pedigree, sires, and dams, between others), which can explain differences in some of the genetic parameters. In general, according to Zaalberg et al. (2023) based on a range of literature reports (Roehe, 1999; Kaufmann et al., 2000; Grandinson et al., 2002), the range of values for the direct and maternal genetic heritability for BW in conventional maternal pig breeds (Landrace and Yorkshire) is between 0.04 to 0.15 and 0.14 to 0.28, respectively. The negative correlation between the direct additive and maternal genetic effects observed was also reported in previous studies (Zhang et al., 2000; Grandinson et al., 2002; Arango et al., 2006; Muns et al., 2013; Alves et al., 2018). This can be explained by the fact that the data structure might not be optimal to properly estimate the correlation between the direct additive and maternal genetic effects or because non-negligible cross-fostering exists (Alves et al., 2018). Cross-fostering information was not known in this study.

Birth weight is an important trait in the swine industry because it is related to piglet survival, behavior, and weight gain (Ayuso et al., 2020; Knol et al., 2022). Low BW is associated with poor gastrointestinal development, which impacts nutrient digestion, absorption, and transportation, implying future health, welfare, and market problems in pigs (Ayuso et al., 2020). Considering the proportion of phenotypic variance explained by the epigenetic effect in the piglet’s BW (0.33 – Table 3), attention should be given to the sow environment during pregnancy, such as weather, sanitary management, and nutrition, as these components can induce epigenetic changes in the progeny (Ibeagha-Awemu and Yu, 2021).

The heritability estimates for the direct genetic effect of the WW in pigs range from 0.01 to 0.22 (Tomiyama et al., 2010; Muns et al., 2013; Dufrasne et al., 2014; Jiao et al., 2014; Alves et al., 2018; Zaalberg et al., 2023) and for the maternal genetic effect from 0.06 to 0.24 (Grandinson et al., 2002; Tomiyama et al., 2010; Jiao et al., 2014; Alves et al., 2018). Similar to previous studies listed above, the heritability of the maternal genetic component reduced from birth to weaning, which reflects the decrease in contribution of the sow effect to the individual body weight of her offspring. The individual’s own genetics is expected to be more important as the pigs age (Alves et al., 2018; Muns et al., 2013). Differences in the heritability estimates for WW among studies, beyond the previous causes mentioned before (i.e., data structure and population), could also be due to the model used. The correlation between the direct additive and maternal genetic effects was not significant in this study for WW (Table 4), as reported in other studies (Alves et al., 2018; Kaufmann et al., 2000; Grandinson et al., 2002). Data structure and cross-fostering can also be reasons for the lack of significance in the correlation between the direct additive and maternal genetic effects for WW.

The weaning process is a stressful situation for both piglets and dams, marked by the time when the piglet leaves the maternal to a different environment and changes the diet, characterizing social and dietary stress (Corbett et al., 2013; Ibeagha-Awemu and Yu, 2021). The WW can be used to assess the success of the lactation phase, and it is the start of the period when the pigs will express their genetic potential without their mother’s direct contributions. Heavy pigs at weaning tend to also be heavy at birth and slaughter age and have faster growth rates, achieving the age to slaughter earlier than lighter pigs (Wolter and Ellis, 2001). Transgenerational epigenetic heritability for WW in Landrace pigs explained more than 30% of the phenotypic variation, and more research is needed to understand the impact of epigenetics due to production management practices. This way, strategies to minimize stress at weaning are important for reducing the potential negative impacts on the pig’s performance.

Beyond growth, desirable body composition characteristics are important breeding goals in the pig industry. The genetic parameters for BF in this study are within the range reported in the literature for body composition traits (0.33–0.65), which are known to have moderate to high direct narrow-sense heritability (Bryner et al., 1992; Singh et al., 2001; Kim et al., 2004; Akanno et al., 2013; Khanal et al., 2019; Willson et al., 2020). Maternal genetic heritability was also a significant variance component in a genetic model for BF in pigs in a previous study (Bryner et al., 1992), with the estimate higher than the one observed in this study (0.11). Bryner et al. (1992) studied Yorkshire pigs and found a negative genetic correlation between direct and maternal genetic components for BF, which was stronger than the estimate observed in this study (−0.51). Negative correlations between direct and maternal genetic effects for BF can happen due to negative genetic association between these effects or data structure and modeling (Bryner et al. ,1992). These results suggest that BF can have a high response to direct selection and that maternal effects should be considered in the genetic models of the studied population.

Beyond the high narrow-sense heritability, BF can be an indicator of the overall body fat and be easily measured through ultrasound, being an important criterion to select pigs for lean carcass (Gozalo-Marcilla et al., 2021). The reduction in carcass fatness can imply higher growth efficiency and lean meat content (Lonergan et al., 2001). In addition, BF is important for the reproductive performance of sows (Roongsitthichai and Tummaruk, 2014). Roongsitthichai and Tummaruk (2014) reviewed the importance of BF for reproductive traits in pigs and reported that gilts with higher BF attained puberty earlier, delivered one more pig, produced heavier piglets with higher growth rates, and had shorter weaning to service intervals. Interestingly, the transgenerational epigenetic heritability for BF also explained more than 30% of the phenotypic variation, similar to BW and WW. Considering the importance of BF for the pig industry and the percentage of phenotypic variance explained by the transgenerational epigenetic effects, more research on this topic is needed.

The variance components and genetic parameters for TNB and NBA were similar in this research and, therefore, they were discussed together. This similarity in the genetic parameters for TNB and NBA was reported before (Akanno et al., 2013; Putz et al., 2015; Ogawa et al., 2019; Paixao et al., 2019; Yang et al., 2023) with both traits being highly genetically correlated (Alves et al., 2018). Reproductive traits usually present low heritability and the range of heritability estimates for TNB and NBA is from 0.01 up to 0.20 (Akanno et al., 2013; Putz et al., 2015; Alves et al., 2018; Ogawa et al., 2019; Paixao et al., 2019; Yang et al., 2023), with most of the estimates varying from 0.06 to 0.12. The repeatability of the TNB and NBA estimates in the literature range from 0.03 to 0.24 (Alves et al., 2018; Ogawa et al., 2019; Paixao et al., 2019; Yang et al., 2023). The heritability and repeatability estimates for TNB and NBA in this study are close to the estimates from previous studies. As previously discussed, the population and data structure are the main reasons for the differences observed in the literature, as the models for these traits tend to be similar.

As expected, and different from what was observed for growth and body composition traits, the transgenerational epigenetic variance and heritability were low for TNB and NBA. We expected that traits measured late in life would be less impacted by epigenetic effects because: 1) the direct genetic effects tend to increase in importance as the animals age; and 2) temporary environmental effects (environmental effects not fully captured by data structure or model) could be the main source of environmental variation. At this point, it is important to remember that transgenerational epigenetic effects are due to epigenetic markers inherited, and they should affect the animal’s life permanently (Heard and Martienssen, 2014). In this context, the source of transgenerational epigenetic variance in the repeatability models should also be considered. Different from single-measurement models in this study, as for BW, WW, and BF, the transgenerational epigenetic variance in the repeatability models for the TNB and NBA were all extracted from the permanent environment effect. The permanent environment accounts for non-additive genetic effects (i.e., dominance, epistasis), as shown by Vitezica et al. (2018), and any other effect with a long-lasting impact on the animal’s life. With this result, it is possible to conclude that the transgenerational epigenetic effect could be treated as a non-additive genetic effect in repeatability models and being accounted for in the genetic models through multiple measurements. This way, it becomes easy to see that including the transgenerational epigenetic effects in the genetic models is more important for traits that can be measured only once in the animal’s life, e.g., BW, WW, and age at first farrowing.

Beyond the repeatability models, another alternative to fit permanent environment effects in a genetic model is by using random regression models (RRM), which also allows for modeling the phenotype (and its components) trajectory through covariance functions (Schaeffer, 2004). In this context, we also expect that RRM can account for the transgenerational epigenetic effects in the permanent environment effects, but no precedents were found in the literature. We recommend such studies with better data structures (e.g., multiple measurements for growth traits early in life) to test this hypothesis in the future. However, the confounding between the transgenerational epigenetic and permanent environment effects and the impacts of this in convergence should be considered, especially for complex models like RRM.

Several differentially methylated regions (DMRs) were reported in the literature for BW (Ayuso et al., 2020), WW (Corbett et al., 2013), BF (Li et al., 2012; Zhang et al., 2016; Wang K. et al., 2022), and TNB and NBA (Hwang et al., 2017; An et al., 2019) in pigs. The DMRs are the most studied epigenetic marks and affect gene expression (Ayuso et al., 2020; Corbett et al., 2013; Li et al., 2012; Zhang et al., 2016; Wang K. et al., 2022; Hwang et al., 2017; An et al., 2019), consequently, can affect the phenotypic variation and be the cause of the epigenetic heritability estimated for the traits investigated in this study. In previous reports, the transgenerational epigenetic heritability was 0.04 for BW in the Pirenaica beef cattle (Varona et al., 2015), and 0.00 to 0.10 for body weight traits in meat quails (Paiva et al., 2018b), and 0.00 to 0.07 in egg traits in meat quails (Paiva et al., 2018a). No previous reports on pigs were found. The species, models, and production systems, among other factors previously mentioned, may be the cause of the difference between the transgenerational epigenetic heritability estimates in the literature and in this study. A reason for different epigenetic heritability across traits could be because epigenetic marks can be linked to QTL with different patterns of linkage disequilibrium (LD) with the traits evaluated. In other words, the epigenetic marks can overlap with or be nearby QTL that are in more LD with a specific trait than others. Therefore, epigenetic heritability is expected to be higher for traits that have more QTL affected by epigenetic marks. One way to find insights into this subject would be by crossing DNA methylation maps with QTL maps for traits with known epigenetic heritability, which was not possible in this study because no whole-genome DNA methylation data was available for the animals used. Future work on this topic is recommended. Following Animal QTLdb, which is a public QTL information repository for several farm species, including pigs (Hu et al., 2022), the development of public DNA methylation maps in pigs would also help to integrate information and understand the causes of the transgenerational epigenetic variance and heritability found in this and future studies. DNA methylation maps in pigs are still scarce (Wang and Kadarmideen, 2019).

One of the main objectives of this study was to compare the prediction results from the best epigenetic model against the genetic or genomic model. As there was no stability in the variance components of the Epi-ssGBLUP, shown by the absence of orthogonality in the models, no decorrelations of the relationships were performed, and no predictions were made for Epi-ssGBLUP. The Epi-ssGBLUP may have been affected by an un-modeled covariance among the epigenetic effects and that of the genetic and residual effects, resulting in drastic changes in the variance components. There were no substantial differences between the variance components from the BLUP and ssGBLUP, suggesting that the variance components from the BLUP model are consistent even though the population evaluated is being genomically selected. Using variance components from Epi-BLUP models, which were stable and more reliable, predictions with the Epi-ssGBLUP model were not performed because making predictions in a model using a relationship matrix different from the one used to estimate the variance components is not recommended (Nilforooshan and Ruíz-Flores, 2022). In the end, the BLUP models were used to compare the predictions against the Epi-BLUP models for similar reasons explained in the Variance Component Estimation section.

By using all animals present in the variance components estimation analysis, no differences are expected in the genetic evaluation, as the correlation between the EBVs estimated in genetic and epigenetic models was high (>0.85; Figure 2). The same could be interpreted for the other random effects in the models, except for the maternal solutions in the BF (−0.22). It is important to highlight that the correlations presented in Figure 2 should not be interpreted as prediction accuracy. The solutions used to compute the Pearson correlations included all phenotypes available, and all solutions were used only to understand the association between the solutions. A formal comparison of EBV prediction, e.g., prediction accuracies, for selection candidates is provided by the results of the LR method (discussed later). Even though the association between solutions was high, the changes in the means of the solutions were significant, showing a different aggregate response for the solutions, except for the permanent environment. Observing that the epigenetic models had a better fit than the genetic models (please see the Model comparison section in the results and discussion), the averages of the solutions for the additive effects of the BW, WW, and BF from the genetic models were underestimated, and the solutions are more dispersed. For the case of the additive solutions for the TNB and NBA and maternal genetic solutions for BW, WW, and BF, the averages from the genetic models are overestimated and more dispersed compared to the epigenetic models. In this case, including the transgenerational epigenetic effects in the genetic models could help properly estimate the average and distribution of the solutions of the mixed model equations, provided that the epigenetic models are better than the genetic models.

The prediction (population) accuracy in pigs usually ranges between 0.05 and 0.82 (reproductive and body composition traits, respectively) (Forni et al., 2011; Hidalgo et al., 2015a; Hidalgo et al., 2015b; Hidalgo et al., 2016; Badke et al., 2014; Fangmann et al., 2017; Oh et al., 2017; Tan et al., 2017; Lee et al., 2020; Zhao et al., 2019; Ardestani et al., 2021; Kjetså et al., 2022; Wang X. et al., 2022; Zhuang et al., 2023). Several factors can affect the prediction accuracies, such as heritability, LD (Meuwissen et al., 2001), population structure, effective population size (Daetwyler et al., 2012), and size of the reference population (Daetwyler et al., 2010). Beyond these factors, the genetic model, accuracy formula used, validation strategy, breed, and the use of commercial crossbreed data also differ between studies and can be a source variation for the prediction accuracy in the pig literature. Using genomic data, Jiao et al. (2014) reported prediction accuracies of 0.194 for BW in Duroc pigs using the BayesA model. Also, using genomic data, Lee et al. (2020) reported prediction accuracies for BW in Yorkshire pigs with the BayesB and BayesC models using panels of different platforms ranging from 0.150 to 0.261. Even though we used only phenotypes and pedigree, the prediction accuracies obtained for the BW were within this range.

A low prediction accuracy (0.098) was reported for WW in Duroc pigs (Jiao et al., 2014), lower than the estimate obtained in this study (∼0.22). The previous authors used the formula from Legarra et al. (2008) to calculate the prediction accuracy and an imputed 60K from a low-density marker panel to perform the genomic predictions. The accuracy formula from Legarra et al. (2008) depends on the correlation between the observed and expected performance divided by the square root of the heritability of the trait, which can be reasonable if there is no selection, and the fixed effects and variance components are properly estimated (Legarra and Reverter, 2018). If there is selection, the accuracy measured by the method proposed by Legarra et al. (2008) underestimates the prediction accuracy, and problems may also arise if the fixed effects and/or variance components are not well estimated (Legarra and Reverter, 2018). In addition, prediction accuracy when using imputed data is also a function of the imputation accuracy, with higher accuracies expected if the imputation accuracy is high (Bolormaa et al., 2015).

Prediction accuracies ranging from 0.285 to 0.527 were observed for BF in Canadian Landrace by using parent average (PA) and several genomic models (Ardestani et al., 2021), in which the lowest prediction accuracy was obtained by PA and the highest by the ssGBLUP model. Ardestani et al. (2021) measured the prediction accuracy as the correlation between (G)EBV and deregressed EBVs. Despite the differences in the accuracy formula and models between our study and Ardestani et al. (2021), the ranges of prediction accuracies were similar. The prediction accuracy for the EBVs of BF in our study was lower than the ones obtained by the genomic models used by Ardestani et al. (2021), but higher than the PA. Jiao et al. (2014) reported a prediction accuracy lower than the one obtained in our study (0.365) using genomic models.

Most studies investigating prediction accuracy in pigs were found in reproductive traits, especially TNB. For NBA, Fangmann et al. (2017) found prediction accuracies ranging from 0.08 to 0.52 using different validation scenarios and accuracy measures in Landrace pigs. In Chinese Yorkshire pigs, Song et al. (2019) reported accuracies ranging from −0.008 to 0.668 for NBA using different SNP panels and validation scenarios. Wang X. et al. (2022) observed prediction accuracies ranging between 0.207 and 0.328 for NBA also in Chinese Yorkshire but based on several machine learning models and genomic BLUP and ssGBLUP. However, it should be noted that the machine learning models may not result in only the predictions of breeding values. The range of accuracy estimates for TNB reported in the literature range from −0.001 to 0.790 (Forni et al., 2011; Hidalgo et al., 2015a; Hidalgo et al., 2015b; Hidalgo et al., 2016; Oh et al., 2017; Song et al., 2019; Kjetså et al., 2022; Wang X. et al., 2022), based on different models, validation strategies, breeds, and other factors as previously mentioned.

In general, prediction accuracy results in pig literature are reported mainly for genomic models, with some studies presenting pedigree-based prediction accuracies and formal comparisons (Forni et al., 2011; Christensen et al., 2012; Tusell et al., 2013; Hidalgo et al., 2015b; Oh et al., 2017; Ardestani et al., 2021). Including genomics in the prediction models tends to increase prediction accuracies because it allows to capture more information to be used to differentiate animals and correct expected relationships (Mendelian sampling) (Georges et al., 2019). Similar or higher prediction accuracies were observed by some authors using only pedigree when compared to genomic relationships (Hidalgo et al., 2015b; Araujo et al., 2023). This may happen when the heritability of the trait is high, pedigree is relatively deep and accurate, there are rare variants with big effects or insufficient genomic information.

Prediction bias calculated as the difference between the (G)EBV estimated in the whole and partial data (Legarra and Reverter, 2018) are scarce in the pig literature. Jang et al. (2023) observed the prediction bias ranging from −0.003 to 0.31 for BF in three terminal pig lines based on different scenarios using preselected variants from whole genome sequence data. The lowest value observed in the BF prediction bias by these authors was the same as the one observed in this study for both BLUP and Epi-BLUP models. Even though only pedigree and phenotypes were used to investigate prediction results in this study, the prediction bias obtained can be considered low, as it ranged between −0.011 and 0.08 (BLUP models for BW and WW, respectively). In the case of the prediction dispersion, which sometimes is called bias, it has been reported as the slope of the regression coefficient of the adjusted phenotype/EBV on the predicted EBV, different from this research that used a deviation from 1. The prediction dispersion (measured as the slope of the regression) in the pig literature ranges between −1.14 and 5.00 (Fangmann et al., 2017; Hidalgo et al., 2015b; Hidalgo et al., 2016; Wang X. et al., 2022), in which the factors that justify the changes in the literature results are similar to those previously stated for the prediction accuracy. The expected/desired value for the dispersion is 1 if measured as the slope of the regression coefficient of the adjusted phenotype/EBV on the predicted EBV or 0 if measured as the deviation from 1 (as performed in this research). Based in the literature cited above, TNB is the trait that presents worse results for dispersion, and it was the trait with the second worse results for dispersion in our study.

To the best of our knowledge, no prediction results for the measures investigated in this research (accuracy, bias, and dispersion of EBV) in genetic models incorporating epigenetic effects have been reported in the literature. This study provides such results for the first time. The similarity between the prediction accuracies, bias, and dispersion in the BLUP and Epi-BLUP models is expected as the transgenerational epigenetic in the Epi-BLUP model is fitted independent of the other random effects, and the orthogonality was present. This shows that including the epigenetic effects in the genetic model does not change the breeding value prediction of focal animals, which suggests that including epigenetic effects in the model is not required from a practical point of view in the current framework (e.g., not including genomic or epigenomic information, not accounting potential interactions). However, it should be noted that the prediction results in this research were not meant to perform model comparison/selection, as this was the purpose of the Model Comparison section, in which the DIC was used as a criterion and showed that epigenetic models have the best fit for all traits. Prediction results were included to show potential advantages in breeding values predicted by epigenetic models over regular animal models, which was not the case probably because of the reasons explained above. Nevertheless, more research is needed on this topic. In addition, using epigenetic genetic effects may still be important for performance/phenotype prediction as it contributes to phenotypic variation.

There is significant transgenerational epigenetic variance for growth, body composition, and reproductive traits in Landrace pigs. Theoretical and analytical work should be done to propose methods for deriving a covariance between epigenetic and additive and residual effects, especially for genomic models where there was non-orthogonality of the effects included. The estimated percentage of epigenetic marks being inherited in Landrace pigs ranged from 10% to 20%, while the reset coefficient was between 0.80 and 0.90, considering the traits evaluated. Transgenerational epigenetic heritability (fraction of the phenotypic variance explained by the transgenerational epigenetic effects) is high for growth and body composition traits (>0.30) and low for reproductive traits (∼0.04) in Landrace pigs. The permanent environment effect can capture the transgenerational epigenetic variances, classifying the epigenetic effect as another non-additive genetic effect affecting phenotypic variation and indicating that it is more important to fit epigenetic effects in single-record traits. EBVs from genetic and epigenetic models are highly correlated at the population level, while statistical differences between the mean values exist. Including epigenetic effects in genetic models does not impact the prediction accuracy, bias, and dispersion of the EBVs of focal individuals (young non-phenotyped selection candidates). This work provides a comprehensive investigation of the impact of epigenetic effects on genetic and genomic models from a quantitative point of view, providing a solid base for future studies in this area.