Multifactorial methods integrating haplotype and epistasis effects for genomic estimation and prediction of quantitative traits

The rapid growth in genomic selection data provides unprecedented opportunities to discover and utilize complex genetic effects for improving phenotypes, but the methodology is lacking. Epistasis effects are interaction effects, and haplotype effects may contain local high-order epistasis effects. Multifactorial methods with SNP, haplotype, and epistasis effects up to the third-order are developed to investigate the contributions of global low-order and local high-order epistasis effects to the phenotypic variance and the accuracy of genomic prediction of quantitative traits. These methods include genomic best linear unbiased prediction (GBLUP) with associated reliability for individuals with and without phenotypic observations, including a computationally efficient GBLUP method for large validation populations, and genomic restricted maximum estimation (GREML) of the variance and associated heritability using a combination of EM-REML and AI-REML iterative algorithms. These methods were developed for two models, Model-I with 10 effect types and Model-II with 13 effect types, including intra- and inter-chromosome pairwise epistasis effects that replace the pairwise epistasis effects of Model-I. GREML heritability estimate and GBLUP effect estimate for each effect of an effect type are derived, except for third-order epistasis effects. The multifactorial models evaluate each effect type based on the phenotypic values adjusted for the remaining effect types and can use more effect types than separate models of SNP, haplotype, and epistasis effects, providing a methodology capability to evaluate the contributions of complex genetic effects to the phenotypic variance and prediction accuracy and to discover and utilize complex genetic effects for improving the phenotypes of quantitative traits.

Introduction

Genomic estimation of variance components and associated heritabilities and genomic prediction for quantitative traits using single nucleotide polymorphism (SNP) markers and mixed models have become a widely used approach for genetic improvement in livestock and crop species. The rapid growth in genomic selection data provides unprecedented opportunities to discover and utilize complex genetic mechanisms, but methodology and computing tools are lacking for investigating complex genetic mechanisms using the approach of genomic estimation and prediction. The integration of global low-order epistasis effects and local high-order epistasis effects contained in haplotypes for genomic estimation and prediction is a step forward for the discovery and application of complex genetic mechanisms to improve the phenotypes of quantitative traits.

The theory of genetic partition of two-locus genotypic values defines four types of epistasis values: additive × additive (A×A), additive × dominance (A×D), dominance × additive (D×A), and dominance × dominance (D×D) epistasis values (Cockerham, 1954; Kempthorne, 1954). The Cockerham method defines each epistasis coefficient as the product of the coefficients of the two interacting effects that each can be an additive or dominance effect (Cockerham, 1954). This definition of epistasis coefficient is the basis for defining epistasis model matrices in terms of the model matrices of additive and dominance effects. Cockerham also defines a pedigree epistasis relationship as the product between the pedigree additive and dominance relationships (Cockerham, 1954), and this definition is the theoretical basis for Henderson’s approach to express epistasis relationship matrices as the Hadamard products of the additive and dominance relationship matrices (Henderson, 1985). Henderson’s Hadamard products for epistasis relationship matrices were suggested for genomic prediction using epistasis effects by replacing the pedigree additive and dominance relationship matrices with the genomic additive and dominance relationship matrices calculated from SNP markers (Su et al., 2012; Muñoz et al., 2014; Vitezica et al., 2017). This genomic version of Henderson’s Hadamard products avoids the use of large epistasis model matrices that can be difficult or impossible to compute but contains intra-locus epistasis effects that are not present in the epistasis model (Martini et al., 2020). For this reason, the genomic version of Henderson’s Hadamard products could be described as approximate genomic epistasis relationship matrices (AGERM). Formulations have been developed to obtain the exact genomic epistasis relationship matrices (EGERM) that remove the intra-locus epistasis effects in AGERM by modifying Henderson’s Hadamard products without creating the epistasis model matrices (Jiang and Reif, 2015; Martini et al., 2016; Jiang and Reif, 2020; Martini et al., 2020). The difference between AGERM and EGERM tends to diminish as the number of SNPs increases (Jiang and Reif, 2020). Henderson’s Hadamard products and hence AGERM are applicable to any order of epistasis effects, and EGERM also has a general formula for any order of epistasis effects (Jiang and Reif, 2020). However, limited tests showed that fourth-order global epistasis contributed virtually nothing to the phenotypic variance but generated considerable computing difficulty (Liang et al., 2021), raising questions about the value of global epistasis effects beyond the third-order. Methods of genomic estimation and prediction of global epistasis effects up to the third-order should have wide-range applications, given that the number of reported epistasis effects lags far behind the number of single-point effects (Carlborg and Haley, 2004; Phillips, 2008; Ritchie and Van Steen, 2018) even though epistasis effects are important genetic effects (Cordell, 2002; Segre et al., 2005; Mackay, 2014). In contrast to the computing difficulty and uncertain impact of global high-order epistasis effects beyond the third-order, local high-order epistasis effects in haplotypes with potentially many SNPs were responsible for the increased accuracy of predicting phenotypic values of certain traits (Liang et al., 2020; Bian et al., 2021). The integration of haplotype and epistasis effects provides an approach to investigate the contributions of global low-order epistasis effects and local high-order epistasis effects to the phenotypic variance and the accuracy of genomic prediction under the same model.

An epistasis GWAS in Holstein cattle showed that intra- and inter-chromosome epistasis effects affected different traits differently, for example, the daughter pregnancy rate was mostly affected by inter-chromosome epistasis effects, whereas milk production traits were mostly affected by intra-chromosome epistasis effects (Prakapenka et al., 2021), and genomic heritability estimates of intra- and inter-chromosome heritabilities for the daughter pregnancy rate using methods in this article showed that inter-chromosome A×A heritability was much higher than the intra-chromosome A×A heritability (Liang et al., 2022). Therefore, dividing pairwise epistasis effects into intra- and inter-chromosome epistasis effects allows the investigation of the contributions of intra- and inter-chromosome pairwise epistasis effects to the phenotypic variance and prediction accuracy.

The purpose of the multifactorial model in this article is to integrate haplotype effects and epistasis effects up to the third-order for genomic estimation of variance components and associated heritabilities, as well as genomic prediction of genetic and phenotypic values of quantitative traits, to provide a general and flexible methodology framework for genomic prediction and estimation using complex genetic mechanisms and to provide methodology details of the EPIHAP computer package that implements the integration of haplotype and epistasis effects (Liang et al., 2021, 2022). The methodology in this article will facilitate the discovery and utilization of global low-order and local high-order epistasis effects relevant to the phenotypic variances and prediction accuracies of quantitative traits, and obtain new knowledge of complex genetic mechanisms underlying quantitative traits.

Materials and methods

Quantitative genetics model with single nucleotide polymorphism, haplotype, and epistasis effects and values

The mixed model with single-SNP additive and dominance effects, haplotype additive effects, and pairwise SNP epistasis effects in this article is based on the quantitative genetics (QG) model resulting from the genetic partition of single-SNP genotypic values (Da et al., 2014; Wang and Da, 2014), haplotype genotypic values (Da, 2015), and pairwise genotypic values (Cockerham, 1954). An advantage of this QG model is the readily available quantitative genetics interpretations of SNP additive and dominance effects, values, and variances; haplotype additive effects, values, and variances; epistasis effects, values, and variances; and the corresponding SNP, haplotype, and epistasis heritability estimates. Two QG models are developed: Model-I with 10 effect types, including SNP additive and dominance effects, haplotype additive effects, and epistasis effects up to the third-order; and Model-II with 13 effect types resulting from replacing the pairwise epistasis effects of Model-I with intra- and inter-chromosome epistasis effects. Detailed descriptions of the effects, values, model matrices, the coding of the model matrices, as well as the precise definition of each term in the two QG models, are provided in Supplementary Text S1 and Supplementary Table S1. With these precise definitions of genetic effects, values, and model matrices in the QG models, a concise multifactorial QG model covering both Model-I and Model-II can be established, that is

$\mathbf{g}=\mathrm{\mu }\mathbf{I}+{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{W}}_{\mathrm{i}}{\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}=\mathrm{\mu }\mathbf{I}+{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{u}}_{\mathrm{i}}(1)$

${\mathbf{u}}_{\mathrm{i}}={\mathbf{W}}_{\mathrm{i}}{\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}(2)$

where ${\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}=$ genetic effects of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type from the original QG model based on genetic partition, ${\mathbf{W}}_{\mathrm{i}}=$ model matrix of ${\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}$, ${\mathbf{u}}_{\mathrm{i}}=$ genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type from the original QG model, and f = number of effect types. For Model-I, subscripts $\mathrm{i}=1,\dots ,10$ represent SNP additive (A), SNP dominance (D), haplotype additive, A×A, A×D, D×D, A×A×A, A×A×D, A×D×D, and D×D×D effects sequentially. For Model-II, subscripts $\mathrm{i}=1,\dots ,13$ represent SNP additive, SNP dominance, haplotype additive, intra-chromosome A×A, intra-chromosome A×D, intra-chromosome D×D, inter-chromosome A×A, inter-chromosome A×D, inter-chromosome D×D, A×A×A, A×A×D, A×D×D, and D×D×D effects sequentially. The variance–covariance matrix of the genetic values of Eqs 1 and 2 is

$\mathbf{G}=\mathrm{v}\mathrm{a}\mathrm{r}\left({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{W}}_{\mathrm{i}}{\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}\right)={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathbf{u}}_{\mathrm{i}}\right)={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathrm{i}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2{\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }(3)$

$\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}\right)={\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2\mathbf{I}(4)$

${\mathbf{G}}_{\mathrm{i}}=\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathbf{u}}_{\mathrm{i}}\right)={\mathbf{W}}_{\mathrm{i}}\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}\right){\mathbf{W}}_{\mathrm{i}}^{\prime }={\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2{\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }(5)$

where ${\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^{\mathit{2}}=\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathrm{\tau }}_{\mathrm{i}\mathrm{j}\mathrm{o}}\right)$ genetic variance of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type under the original QG model is common to all individuals (all j values). It is of note that ${\mathbf{W}}_{\mathrm{i}}{{\mathbf{W}}_{\mathrm{i}}}^{\prime }$ is not a genomic relationship matrix but is the primary information for calculating each genomic relationship matrix. The structure of the G matrix of Eqn. 3 assumes independence between the genetic values of different effect types. However, the GBLUP values of different effect types using the G matrix of Eqn. 3 could be correlated. Under the Hardy–Weinberg equilibrium (HWE) and LE assumptions, additive, dominance, and epistasis effects are independent of each other (Cockerham, 1954; Kempthorne, 1954). For genome-wide SNPs, the LE assumption generally does not hold for closely linked loci, and nonzero Hardy–Weinberg disequilibrium (HWD) may exist numerically. These and other unknown factors in real data may result in the existence of correlations between different effect types. Haplotype additive values are correlated with SNP additive effects because a haplotype additive value is the sum of all SNP additive values and an epistasis value within the haplotype block plus a potential haplotype loss (Da et al., 2016). In two recent haplotype studies for genomic prediction, the integration of SNP and haplotype effects increased the prediction accuracy for four of the seven traits in the human study (Liang et al., 2020) and for three of the eight traits in the swine study (Bian et al., 2021), showing that SNP and haplotype additive values compensated each other for prediction accuracy and that the correlation between SNP and haplotype additive values were incomplete for those traits. The correlation between haplotype and epistasis values can be complex. The correlation should be nonexistent if the A×A values are inter-chromosome A×A values or intra-chromosome A×A values involving distal SNPs not covered by the haplotypes, but the correlation could be strong if the A×A values are intra-chromosome A×A values involving proximal SNPs covered by the haplotypes.

The reparametrized and equivalent quantitative genetics model for genomic estimation and prediction

Genomic relationship matrices are used for genomic estimation and prediction, and the use of genomic relationship matrices results in a reparametrized and equivalent model of the original QG model for genetic values, to be referred to as the RE-QG model, where “reparametrized” refers to the reparameterization of the genetic effects, model matrix, and genetic variance of each effect type; and “equivalent” refers to the requirement of the same first and second moments for the original QG model (Eqs 1–5) and the RE-QG model. This RE-QG model of genetic values can be expressed as

$\mathbf{g}=\mathrm{\mu }\mathbf{I}+{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{T}}_{\mathrm{i}}{\mathbf{\tau }}_{\mathrm{i}}=\mathrm{\mu }\mathbf{I}+{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{u}}_{\mathrm{i}}(6)$

$\mathbf{G}=\mathrm{v}\mathrm{a}\mathrm{r}\left({\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{u}}_{\mathrm{i}}\right)={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathrm{i}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{T}}_{\mathrm{i}}{\mathbf{T}}_{\mathrm{i}}^{\prime }={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2{\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }(7)$

where

${\mathbf{\tau }}_{\mathrm{i}}=\sqrt{{\mathrm{k}}_{\mathrm{i}}}{\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}=\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}{\mathrm{i}}^{\mathrm{t}\mathrm{h}}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{e}(8)$

${\mathbf{T}}_{\mathrm{i}}={\mathbf{W}}_{\mathrm{i}}/\sqrt{{\mathrm{k}}_{\mathrm{i}}}=\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{x}\mathrm{o}\mathrm{f}{\mathbf{\tau }}_{\mathrm{i}}(9)$

${\mathrm{\sigma }}_{\mathrm{i}}^{\mathit{2}}=\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathrm{\tau }}_{\mathrm{i}\mathrm{j}}\right)=\mathrm{t}\mathrm{r}\left({\mathbf{G}}_{\mathrm{i}}\right)/\mathrm{n}={\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}/\mathrm{n}={\mathrm{k}}_{\mathrm{i}}{\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2(10)$

= variance of the genetic effects of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type common to all individuals= average variance of all individuals for the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type

${\mathbf{u}}_{\mathrm{i}}={\mathbf{T}}_{\mathrm{i}}{\mathbf{\tau }}_{\mathrm{i}}={\mathbf{W}}_{\mathrm{i}}{\mathbf{\tau }}_{\mathrm{i}\mathrm{o}}=\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{v}\mathrm{a}\mathrm{l}\mathrm{u}\mathrm{e}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}{\mathrm{i}}^{\mathrm{t}\mathrm{h}}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{e}(11)$

${\mathbf{G}}_{\mathrm{i}}=\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathbf{u}}_{\mathrm{i}}\right)={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{T}}_{\mathrm{i}}{\mathbf{T}}_{\mathrm{i}}^{\prime }={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}}={\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2{\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }(12)$

= variance–covariance matrix of the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type

${\mathbf{S}}_{\mathrm{i}}={\mathbf{T}}_{\mathrm{i}}{\mathbf{T}}_{\mathrm{i}}^{\prime }={\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }/{\mathrm{k}}_{\mathrm{i}}=\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{o}\mathrm{m}\mathrm{i}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{h}\mathrm{i}\mathrm{p}\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{x}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}{\mathrm{i}}^{\mathrm{t}\mathrm{h}}\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{e}(13)$

${\mathrm{k}}_{\mathrm{i}}=\mathrm{t}\mathrm{r}\left({\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }\right)/\mathrm{n}=\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\mathrm{o}\mathrm{n}\mathrm{a}\mathrm{l}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{o}\mathrm{f}{\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }.(14)$

Equations 8–10 are the reparametrization of the genetic effects, model matrices, and genetic variances of the original QG model, whereas Eqs 11 and 12 show the genetic values and the variance–covariance matrix of the genetic values are the same under the RE-QG and QG models. In Eq.10, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}$ = the genetic variance of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type = the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ diagonal element of the ${\mathbf{G}}_{\mathrm{i}}$ matrix defined by Eq. 12. The ${\mathrm{k}}_{\mathrm{i}}$ formula of Eq. 14 as the average of the diagonal elements of ${\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }$ was originally proposed for genomic additive relationships (Hayes and Goddard, 2010) and was used for genomic dominance relationships (Da et al., 2014; Wang and Da, 2014), haplotype additive genomic relationships (Da, 2015), and pairwise epistasis genomic relationships (Vitezica et al., 2017). The need for this RE-QG model is due to the use of the genomic relationship matrices (e.g., Eq. 13) because the QG model does not contain genomic relationship matrices (Eq. 3). Detailed notations of the QG model of Eqs 1–5 in reference to the RE-QG model described by Eqs 6–14 are summarized in Supplementary Table S1.

The formula of the genomic relationship matrix (${\mathbf{S}}_{\mathrm{i}}$ of Eq. 13) is based on the model matrix of each effect type and can be difficult or impossible to compute if epistasis model matrices are used. This computing difficulty of epistasis model matrices is removed by calculating ${\mathbf{S}}_{\mathrm{i}}$ based on the model matrices of SNP additive and dominance effects without creating the epistasis model matrices using either AGERM or EGERM. AGERM refers to the genomic version of Henderson’s Hadamard products between pedigree additive and dominance relationship matrices (Henderson, 1985), with the pedigree additive and dominance relationship matrices replaced by the genomic additive and dominance relationship matrices (Su et al., 2012; Muñoz et al., 2014; Vitezica et al., 2017). AGERM contains intra-locus epistasis that should not exist (Martini et al., 2020), and EGERM removes intra-locus epistasis from AGERM based on products between genomic additive and dominance relationship matrices (Jiang and Reif, 2020; Martini et al., 2020).

The QG and RE-QG models have the same prediction accuracy due to the equivalence between these two models. The genetic values (${\mathbf{u}}_{\mathrm{i}}$, Eqs 2, 11) and the variance–covariance matrix of the genetic values (${\mathbf{G}}_{\mathrm{i}}$, Equations 5 and 12) under the QG and RE-QG models are identical, although ${\mathbf{u}}_{\mathrm{i}}$ and ${\mathbf{G}}_{\mathrm{i}}$ have different expressions under the QG and RE-QG models. Consequently, the QG model without using genomic relationship matrices and the RE-QG model using genomic relationship matrices have identical accuracy of genomic prediction. The choice of the ${\mathrm{k}}_{\mathrm{i}}$ formula for defining the genomic relationship matrix does not affect the accuracy of genomic prediction but affects the interpretation and application of the genetic variance and genomic relationships for each effect type. Since the interpretation of each genetic variance is a focus, whereas the interpretation of the genomic relationships is not a focus in this study, the interpretation of the genetic variance and associated heritability is the consideration in choosing the ${\mathrm{k}}_{\mathrm{i}}$ formula of Eq.14.

The RE-QG model using genomic relationships (Equations 6–14) has two major advantages over the QG model without using genomic relationship matrices (Equations 1–5), although the two models have the same prediction accuracy. First, the use of genomic relationships, originally proposed for genomic additive relationships (VanRaden, 2008), provides a genomic version of the traditional theory and methods of best linear unbiased prediction (BLUP) that uses pedigree relationships, and this genomic version can utilize a wealth of BLUP-based theory, methods, and computing strategies. Second, the genetic variance of the genetic effects of each effect type under the RE-QG model can be used for estimating genomic heritability, whereas the genetic variance of the genetic effects under the QG model cannot be used for estimating genomic heritability. With the ${\mathrm{k}}_{\mathrm{i}}$ value defined by Eq. 14, the variance of the genetic effects of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type, ${\mathrm{\sigma }}_{\mathrm{i}}^2={\mathrm{k}}_{\mathrm{i}}{\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2$ (Eq.10), has the unique interpretation as the average variance of the genotypic values of all individuals and is a common variance to all individuals. Moreover, ${\mathrm{\sigma }}_{\mathrm{i}}^2={\mathrm{k}}_{\mathrm{i}}{\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2$ is unaffected by the number of levels for each effect type, unless the number of levels such as the number of SNPs is too small to provide sufficient coverage of the genome (Da et al., 2014; Tan et al., 2017; Liang et al., 2020). In contrast, the QG model does not have a method to estimate genetic variance components for calculating genomic heritabilities because ${\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2$ is an inverse function of the number of effect levels. As the number of effect levels such as the number of SNPs increases or decreases, the value of each element in ${\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }$ changes in the same direction and the ${\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2$ estimate changes in the opposite direction, that is, as the number of effect levels increases or decreases, ${\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2$ decreases or increases. Consequently, the ${\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2$ estimate does not have a unique interpretation and cannot be used for estimating genomic heritability (Da et al., 2014). Moreover, the variance of the genetic value of an individual $\left({\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2{\left({\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }\right)}^{\mathrm{j}\mathrm{j}}\right)$ cannot be used for calculating genomic heritability because of the individual specificity of the ${\left({\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }\right)}^{\mathrm{j}\mathrm{j}}$ values, as shown as follows.

The exact relationship between the genetic variance for the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual under the RE-QG model and the QG model can be described based on the ${\mathbf{G}}_{\mathrm{i}}$ matrix defined by Eq. 12:

${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}=\mathrm{V}\mathrm{a}\mathrm{r}\left({\mathrm{u}}_{\mathrm{i}\mathrm{j}}\right)={\mathrm{\sigma }}_{\mathrm{i}}^2{\left({\mathbf{S}}_{\mathrm{i}}\right)}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{i}\mathrm{o}}^2{\left({\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }\right)}^{\mathrm{j}\mathrm{j}}(15)$

where ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}=$ the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ diagonal element of the ${\mathbf{G}}_{\mathrm{i}}$ matrix defined by Eq.12 = the genetic variance of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for the genotypic value of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type, and ${\mathrm{u}}_{\mathrm{i}\mathrm{j}}$ = the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ element of ${\mathbf{u}}_{\mathrm{i}}$ defined by Eq.11. Equation 15 shows that different individuals do not have a common variance of the genetic values $\left({\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}\right)$ unless all diagonal elements of ${\mathbf{S}}_{\mathrm{i}}$ or ${\mathbf{W}}_{\mathrm{i}}{\mathbf{W}}_{\mathrm{i}}^{\prime }$ are identical, which could not happen with genome-wide SNP data in the absence of identical twins because genome-wide SNPs have a high degree of individual specificity. Consequently, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}$ is not a common variance to all individuals and cannot be used for calculating the genomic heritability of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type. In contrast, ${\mathrm{\sigma }}_{\mathrm{i}}^2$ of Eq.10 under the RE-QG model as the average variance of the genotypic values of all individuals is common to all individuals and can be used for calculating the heritability of each effect type. For the example of Model-I, the exact genetic interpretation of ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}$ is ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{a}\mathrm{j}}^2=$ the variance of the genomic additive (breeding) value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=1$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{d}\mathrm{j}}^2=$ the variance of the genomic dominance value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=2$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{a}\mathrm{h}\mathrm{j}}^2=$ the variance of the genomic haplotype additive value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=3$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{a}\mathrm{a}\mathrm{j}}^2=$ the variance of the A×A value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=4$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{a}\mathrm{d}\mathrm{j}}^2=$ the variance of the A×D value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=5$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{d}\mathrm{d}\mathrm{j}}^2=$ the variance of the D×D value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=6$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{a}\mathrm{a}\mathrm{a}\mathrm{j}}^2=$ the variance of the A×A×A value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=7$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{a}\mathrm{a}\mathrm{d}\mathrm{j}}^2=$ the variance of the A×A×D value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=8$, ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{a}\mathrm{d}\mathrm{d}\mathrm{j}}^2=$ the variance of the A×D×D value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=9$, and ${\mathrm{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\mathrm{\sigma }}_{\mathrm{d}\mathrm{d}\mathrm{d}\mathrm{j}}^2=$ the variance of the D×D×D value of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual for $\mathrm{i}=10$. These genetic interpretations, along with those for intra- and inter-chromosome pairwise epistasis effects of Model-II under the QG and RE-QG models, are summarized in Supplementary Table S1.

Results and discussion

The multifactorial model of phenotypic values

Based on the RE-QG model of Eqs 6–14, the multifactorial model for phenotypic values is

$\begin{array}{l}\mathbf{y}=\mathbf{X}\mathbf{b}+\mathbf{Z}\mathbf{g}+\mathbf{e}=\mathbf{X}\mathbf{b}+\mathbf{Z}{\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{T}}_{\mathrm{i}}{\mathbf{\tau }}_{\mathrm{i}}+\mathbf{e}\\ =\mathbf{X}\mathbf{b}+\mathbf{Z}{\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{u}}_{\mathrm{i}}+\mathbf{e}\end{array}(16)$

$\begin{array}{l}\mathbf{V}=\mathbf{Z}\mathbf{G}{\mathbf{Z}}^{\prime }+{\mathrm{\sigma }}_{\mathrm{e}}^2{\mathbf{I}}_{\mathrm{N}}=\mathbf{Z}\left({\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathrm{i}}\right){\mathbf{Z}}^{\prime }+{\mathrm{\sigma }}_{\mathrm{e}}^2{\mathbf{I}}_{\mathrm{N}}\\ =\mathbf{Z}\left({\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{T}}_{\mathrm{i}}{\mathbf{T}}_{\mathrm{i}}^{\prime }\right){\mathbf{Z}}^{\prime }+{\mathrm{\sigma }}_{\mathrm{e}}^2{\mathbf{I}}_{\mathrm{N}}=\mathbf{Z}\left({\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}}\right){\mathbf{Z}}^{\prime }+{\mathrm{\sigma }}_{\mathrm{e}}^2{\mathbf{I}}_{\mathrm{N}}\end{array}(17)$

where y = N×1 column vector of phenotypic observations, Z = $\mathrm{N}\times \mathrm{n}$ incidence matrix allocating phenotypic observations to each individual = identity matrix for one observation per individual (N = n), N = number of observations, n = number of individuals, $\mathbf{b}=\mathrm{c}\times \mathrm{l}$ column vector of fixed effects such as heard-year-season in dairy cattle, c = number of fixed effects, $\mathbf{X}=\mathrm{N}\times \mathrm{c}$ model matrix, b,e = $\mathrm{N}\times 1$ column vector of random residuals, ${\mathrm{\sigma }}_{\mathrm{e}}^2$ = residual variance, and $\mathbf{G}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathrm{i}}$ (Eq. 7). The phenotypic values $\left(\mathbf{y}\right)$ are assumed to follow a normal distribution with mean $\mathbf{X}\mathbf{b}$ and variance–covariance matrix of V. The methods described below for genomic estimation and prediction are based on the conditional expectation (CE) method, which is more efficient computationally than the methods based on mixed-model equations (MME) when the number of genetic effects is greater than the number of individuals (Da et al., 2014; Da, 2015).

For Model-I with 10 effect types, the genomic epistasis relationship matrices can be calculated using either EGERM or AGERM. However, EGERM or AGERM did not consider intra- and inter-chromosome genomic epistasis relationship matrices that are required by Model-II with 13 effect types. This research derives intra- and inter-chromosome genomic epistasis relationship matrices for both EGERM and AGERM.

Intra- and inter-chromosome genomic epistasis relationship matrices

The main derivation of the intra- and inter-chromosome genomic epistasis relationship matrices is the partition of the numerator of a genomic epistasis relationship matrix into intra- and inter-chromosome numerators. The first step is to derive the intra-chromosome numerator, and the second step is to derive the inter-chromosome numerator as the difference between the whole-genome numerator and the intra-chromosome numerator. The last step is to divide the intra-chromosome numerator by the average of the diagonal elements of the intra-chromosome numerator and to divide the inter-chromosome numerator by the average of the diagonal elements of the inter-chromosome numerator. Using this procedure, intra- and inter-chromosome epistasis relationship matrices were derived for both EGERM and AGERM (Supplementary Text S1).

Genomic best linear unbiased prediction and reliability

Based on the multifactorial genetic model of Eqs 16 and 17, the GBLUP of the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type (${\widehat{\mathbf{u}}}_{\mathrm{i}}$) and the best linear unbiased estimator (BLUE) or generalized least squares (GLS) estimator of fixed effect ($\widehat{\mathbf{b}}$) are

$\begin{array}{cc}{\widehat{\mathbf{u}}}_{\mathrm{i}}={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}}{\mathbf{Z}}^{\prime }{\mathbf{V}}^{-1}\left(\mathbf{y}-\mathbf{X}\widehat{\mathbf{b}}\right)={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}}{\mathbf{Z}}^{\prime }\mathbf{P}\mathbf{y},& \mathrm{i}=1,...,\end{array}\mathrm{f}(18)$

$\widehat{\mathbf{b}}={\left({\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{y}(19)$

where $\mathbf{P}={\mathbf{V}}^{-1}-{\mathbf{V}}^{-1}\mathbf{X}{\left({\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}\mathbf{X}\right)}^{-}{\mathbf{X}}^{\prime }{\mathbf{V}}^{-1}$. The GBLUP of total genetic values of the n individuals is the summation of all types of genetic values:

$\widehat{\mathbf{g}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\widehat{\mathbf{u}}}_{\mathrm{i}}.(20)$

Reliability of GBLUP is the squared correlation between the GBLUP of a type of genetic value and the unobservable true genetic value being predicted by the GBLUP. The expected accuracy of predicting the genetic values by the GBLUP is the square root of reliability or the correlation between the GBLUP of a type of genetic effect and the unobservable true genetic effects being predicted by the GBLUP. In the absence of validation studies for observed prediction accuracy, reliability or the expected prediction accuracy is the measure of prediction accuracy of the GBLUP. The reliability of the GBLUP of the total genetic value (Eq. 2) of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individual is

${\mathrm{R}}_{\mathrm{g}\mathrm{j}}^2={\left[\mathbf{G}\left({\mathbf{Z}}^{\prime }\mathbf{P}\mathbf{Z}\right)\mathbf{G}\right]}^{\mathrm{j}\mathrm{j}}/{\mathbf{G}}^{\mathrm{j}\mathrm{j}}(21)$

where $\mathbf{G}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathbf{i}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{T}}_{\mathrm{i}}{\mathbf{T}}_{\mathrm{i}}^{\prime }={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}}$ (Eq. 7), ${\mathbf{G}}^{\mathrm{j}\mathrm{j}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}}^{\mathrm{j}\mathrm{j}}$, and subscript or superscript jj denotes the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ diagonal element. The reliability formula for any or a combination of genetic values can be readily derived from Eq. 21, for example, the reliability of ${\widehat{\mathbf{u}}}_3$ (GBLUP of haplotype additive values) is obtained from Eq. 21 by deleting all terms except ${\mathbf{G}}_3\left({\mathbf{Z}}^{\prime }\mathbf{P}\mathbf{Z}\right){\mathbf{G}}_3$ in the numerator and ${\mathrm{\sigma }}_3^2{\mathbf{S}}_3^{\mathrm{j}\mathrm{j}}$ in the denominator, with changes in the V and P matrices accordingly.

Calculation of genomic best linear unbiased prediction and reliability for individuals with and without phenotypic observations separately

Two strategies are available for calculating GBLUP and the reliability of Eqs 20 and 21. Strategy-1 is a one-step strategy that includes all individuals with and without phenotypic observations in the same system of equations so that GBLUP and reliability are calculated simultaneously for all individuals. This strategy essentially augments the mixed model for individuals with phenotypic observations with a set of null equations consisting of “0”s but uses each genomic relationship matrix for all individuals, and these null equations and the use of the relationship matrix for all individuals do not affect the GBLUP, reliability, and heritability of individuals with phenotypic observations. The advantage of this one-step strategy is the simplicity of data preparation. For example, for a k-fold cross validation study, the phenotypic input file only needs to have k columns of the trait observations, with one column for each validation where the phenotypic observations for the validation individuals are set as “missing,” and the X and Z model matrices for the “missing” observations are set to zero. With this strategy, the genotypic data need to be processed only once. As the number of traits increases for validation studies, this one-step strategy becomes more appealing due to the savings in data preparation work. This strategy has been implemented in our computing tools of GVCBLUP (Wang et al., 2014), GVCHAP (Prakapenka et al., 2020), and EPIHAP (Liang et al., 2021, 2022). However, when the number of validation individuals or individuals without phenotypic values is large, each genomic relationship matrix (${\mathbf{S}}_{\mathrm{i}}$ matrix) is large, and the one-step strategy becomes more difficult as the number of individuals increases.

For large numbers of individuals without phenotypic observations, calculating GBLUP for individuals with and without phenotypic values separately is more efficient computationally than calculating GBLUP for all individuals in the same system of equations by applying Henderson’s BLUP for animals without phenotypic observations (Henderson, 1977) to GBLUP. Let ${\mathrm{n}}_1$ = number of individuals with phenotypic observations, ${\mathrm{n}}_0$ = number of individuals without phenotypic observations, $\mathrm{n}={\mathrm{n}}_1+{\mathrm{n}}_0$, and let the ${\mathbf{S}}_{\mathrm{i}}$ matrix be partitioned as

$\begin{array}{cc}{\mathbf{S}}_{\mathrm{i}}=\left[\begin{array}{cc}{\mathbf{S}}_{\mathrm{i}11}& {\mathbf{S}}_{\mathrm{i}10}\\ {\mathbf{S}}_{\mathrm{i}01}& {\mathbf{S}}_{\mathrm{i}00}\end{array}\right],& \mathrm{i}=1,\dots ,\mathrm{f}\end{array}(22)$

where ${\mathbf{S}}_{\mathrm{i}11}={\mathrm{n}}_1\times {\mathrm{n}}_1$ genomic relationship matrix of the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type for individuals with phenotypic observations, ${\mathbf{S}}_{\mathrm{i}01}={\mathrm{n}}_0\times {\mathrm{n}}_1=$ genomic relationship matrix of the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type between individuals without phenotypic observations and individuals with phenotypic observations, ${\mathbf{S}}_{\mathrm{i}10}={\mathbf{S}}_{\mathrm{i}01}^{\prime }={\mathrm{n}}_1\times {\mathrm{n}}_0=$ genomic relationship matrix between individuals with phenotypic observations and individuals without phenotypic observations, and ${\mathbf{S}}_{\mathrm{i}00}$ = ${\mathrm{n}}_0\times {\mathrm{n}}_0$ genomic relationship matrix of the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type for individuals without phenotypic observations. In Eqs 16 and 17, $\mathbf{y}={\mathbf{y}}_1$, the Z matrix needs to be changed to $\mathbf{Z}=\left[\begin{array}{cc}{\mathbf{Z}}_1& \mathbf{0}\end{array}\right]$, the ${\mathbf{u}}_{\mathrm{i}}$ vector partitioned as ${\mathbf{u}}_{\mathrm{i}}={\left[{\mathbf{u}}_{\mathrm{i1}}^{\prime },{\mathbf{u}}_{\mathrm{i0}}^{\prime }\right]}^{\prime }$, and the g vector partitioned as $\mathbf{g}=\left[{\mathbf{g}}_1^{\prime },{\mathbf{g}}_0^{\prime }\right]\prime$, where ${\mathbf{Z}}_1=\mathrm{N}\times {\mathrm{n}}_1$ incidence matrix allocating phenotypic observations to individuals with phenotypic observations, $\mathbf{0}=\mathrm{N}\times {\mathrm{n}}_0$ incidence matrix with elements “0” connecting phenotypic observations to individuals without phenotypic observations. With these changes and Eq. 22, the V matrix of Eq. (17) can be re-written as

$\mathbf{V}={\mathbf{Z}}_1\left({\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathrm{i}}\right){\mathbf{Z}}_1^{\prime }+{\mathrm{\sigma }}_{\mathrm{e}}^2{\mathbf{I}}_{\mathrm{N}}={\mathbf{Z}}_1\left({\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}11}\right){\mathbf{Z}}_1^{\prime }+{\mathrm{\sigma }}_{\mathrm{e}}^2{\mathbf{I}}_{\mathrm{N}}(23)$

and the GBLUP and reliability for individuals with and without phenotypic observations can be calculated as

$\begin{array}{cc}{\widehat{\mathbf{u}}}_{\mathrm{i}1}={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}11}{\mathbf{Z}}_1^{\prime }{\mathbf{V}}^{-1}\left({\mathbf{y}}_1-\mathbf{X}\widehat{\mathbf{b}}\right)={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}11}{\mathbf{Z}}_1^{\prime }\mathbf{P}{\mathbf{y}}_1,& \mathrm{i}=1,...,\mathrm{f}\end{array}(24)$

${\widehat{\mathbf{g}}}_1={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\widehat{\mathbf{u}}}_{\mathrm{i}1}(25)$

${\mathrm{R}}_{\mathrm{g}1\mathrm{j}}^2={\left[{\mathbf{G}}_{11}\left({\mathbf{Z}}_1^{\prime }\mathbf{P}{\mathbf{Z}}_1^{\prime }\right){\mathbf{G}}_{11}\right]}^{\mathrm{j}\mathrm{j}}/{\mathbf{G}}_{11}^{\mathrm{j}\mathrm{j}}(26)$

$\begin{array}{cc}{\widehat{\mathbf{u}}}_{\mathrm{i}0}={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}01}{\mathbf{Z}}_1^{\prime }{\mathbf{V}}^{-1}\left({\mathbf{y}}_1-\mathbf{X}\widehat{\mathbf{b}}\right)={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}01}{\mathbf{Z}}_1^{\prime }\mathbf{P}{\mathbf{y}}_1,& \mathrm{i}=1,...,\mathrm{f}\end{array}(27)$

$\begin{array}{c}={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}01}{\mathbf{S}}_{\mathrm{i}11}^{-1}{\mathbf{S}}_{\mathrm{i}11}{\mathbf{Z}}_1^{\prime }\mathbf{P}{\mathbf{y}}_1={\mathbf{G}}_{\mathrm{i}01}{\mathbf{G}}_{\mathrm{i}11}^{-1}{\mathbf{G}}_{\mathrm{i}11}{\mathbf{Z}}_1^{\prime }\mathbf{P}{\mathbf{y}}_1={\mathbf{G}}_{\mathrm{i}01}{\mathbf{G}}_{\mathrm{i}11}^{-1}{\widehat{\mathbf{u}}}_{\mathrm{i}1}\end{array}(28)$

${\widehat{\mathbf{g}}}_0={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\widehat{\mathbf{u}}}_{\mathrm{i}0}(29)$

${\mathrm{R}}_{\mathrm{g}0\mathrm{j}}^2={\left[{\mathbf{G}}_{01}\left({\mathbf{Z}}_1^{\prime }\mathbf{P}{\mathbf{Z}}_1^{\prime }\right){\mathbf{G}}_{10}\right]}^{\mathrm{j}\mathrm{j}}/{\mathbf{G}}_{00}^{\mathrm{j}\mathrm{j}}(30)$

where ${\widehat{\mathbf{u}}}_{\mathrm{i}1}={\mathrm{n}}_1\times 1$ column vector of the GBLUP of the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type for individuals with phenotypic observations, ${\widehat{\mathbf{g}}}_1={\mathrm{n}}_1\times 1$ column vector of the GBLUP of the total genetic values for individuals with phenotypic observations, ${\mathrm{R}}_{\mathrm{g}1\mathrm{j}}^2=$ reliability for the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$individuals with phenotypic observations, ${\widehat{\mathbf{u}}}_{\mathrm{i}0}={\mathrm{n}}_0\times 1$ column vector of the GBLUP of the genetic values of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type for individuals without phenotypic observations, ${\widehat{\mathbf{g}}}_0={\mathrm{n}}_0\times 1$ column vector of the GBLUP of the total genetic values for individuals without phenotypic observations, ${\mathrm{R}}_{\mathrm{g}0\mathrm{j}}^2=$ reliability for the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ individuals without phenotypic observations, ${\mathbf{G}}_{11}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathbf{i}11}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}11}$, ${\mathbf{G}}_{01}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathbf{i}01}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}01}$, ${\mathbf{G}}_{10}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{G}}_{\mathbf{i}10}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{S}}_{\mathrm{i}10}$, ${\mathbf{G}}_{11}^{\mathrm{j}\mathrm{j}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{S}}_{\mathrm{i}11}^{\mathrm{j}\mathrm{j}}{\mathrm{\sigma }}_{\mathrm{i}}^2$, and ${\mathbf{G}}_{00}^{\mathrm{j}\mathrm{j}}={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{S}}_{\mathrm{i}00}^{\mathrm{j}\mathrm{j}}{\mathrm{\sigma }}_{\mathrm{i}}^2$.

Equations 27 and 28 yield identical results if ${\mathbf{S}}_{\mathrm{i}11}^{-1}$ exists. However, when the number of individuals is greater than the number of effect levels, such as the number of SNPs, ${\mathbf{S}}_{\mathrm{i}11}^{-1}$ in Eq. 28 does not exist, and Eq. 27 still can calculate the GBLUP. The usefulness of Eq. 28 is that it shows the GBLUP of individuals without phenotypic observations is the regression of the genetic values of individuals without phenotypic observations on the genetic values of individuals with phenotypic observations. The advantage of Eq. 27 is that it does not calculate ${\mathbf{S}}_{\mathrm{i}11}^{-1}$ and hence is unaffected by the singularity of ${\mathbf{S}}_{\mathrm{i}11}$. Therefore, Eq. 27 is recommended for calculating GBLUP for individuals without phenotypic observation when the number of such individuals is large. The GBLUP calculations of Eqs 24, 27, and 28 do not involve the genomic relationship matrix among individuals without phenotypic observations ${\mathbf{S}}_{\mathrm{i}00}$, which is much larger than ${\mathbf{S}}_{\mathrm{i}11}$ when ${\mathrm{n}}_1$ is much larger than ${\mathrm{n}}_0$. The reliability calculation for individuals without phenotypic observations (Eq. 30) only uses the diagonal elements of ${\mathbf{S}}_{\mathrm{i}00}$ and not the entire ${\mathbf{S}}_{\mathrm{i}00}$.

Advantage of the integrated model over separate models

The multifactorial model of Eqs 16 and 17 integrating SNP, haplotype, and epistasis effects has the advantage of using more effect types and assessing each effect type based on the phenotypic values adjusted for all remaining effect types over separate models for SNP, haplotype, and epistasis effects that do not have a mechanism to adjust for effect types not in the model, and each uses a smaller number of genetic effects in the model.

This advantage of the multifactorial model assessing each effect type based on the phenotypic values adjusted for all remaining effect types can be shown using the MME version of the GBLUP for the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type:

$\begin{array}{l}{\widehat{\mathbf{u}}}_{\mathrm{i}}={\left({\mathbf{Z}}_{\mathrm{i}}^{\prime }{\mathbf{Z}}_{\mathrm{i}}+{\mathbf{G}}_{\mathrm{i}}^{-1}\right)}^{-1}\left[{\mathbf{Z}}_{\mathrm{i}}^{\prime }\mathbf{y}-\left({\mathbf{Z}}_{\mathrm{i}}^{\prime }\mathbf{X}\widehat{\mathbf{b}}+{\sum }_{\begin{array}{l}\mathrm{j}=1\\ \mathrm{j}\ne \mathrm{i}\end{array}}^{\mathrm{f}}{\mathbf{Z}}_{\mathrm{i}}^{\prime }{\mathbf{Z}}_{\mathrm{j}}{\widehat{\mathbf{u}}}_{\mathrm{j}}\right)\right]\\ ={\left({\mathbf{Z}}_{\mathrm{i}}^{\prime }{\mathbf{Z}}_{\mathrm{i}}+{\mathbf{G}}_{\mathrm{i}}^{-1}\right)}^{-1}{\mathbf{Z}}_{\mathrm{i}}^{\prime }\left(\mathbf{y}-\mathbf{X}\widehat{\mathbf{b}}-{\sum }_{\begin{array}{l}\mathrm{j}=1\\ \mathrm{j}\ne \mathrm{i}\end{array}}^{\mathrm{f}}{\mathbf{Z}}_{\mathrm{j}}{\widehat{\mathbf{u}}}_{\mathrm{j}}\right)={\left({\mathbf{Z}}_{\mathrm{i}}^{\prime }{\mathbf{Z}}_{\mathrm{i}}+{\mathbf{G}}_{\mathrm{i}}^{-1}\right)}^{-1}{\mathbf{Z}}_{\mathrm{i}}^{\prime }{\mathbf{y}}_{\mathrm{b}\mathrm{u}}^{*}\end{array}(31)$

$\begin{array}{l}\widehat{\mathbf{b}}={\left({\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-}\left({\mathbf{X}}^{\prime }\mathbf{y}-{\mathbf{X}}^{\prime }{\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{Z}}_{\mathrm{i}}{\widehat{\mathbf{u}}}_{\mathrm{i}}\right)\\ ={\left({\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-}{\mathbf{X}}^{\prime }\left(\mathbf{y}-{\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{Z}}_{\mathrm{i}}{\widehat{\mathbf{u}}}_{\mathrm{i}}\right)={\left({\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-}{\mathbf{X}}^{\prime }{\mathbf{y}}_{\mathrm{u}}^{*}\end{array}(32)$

where ${\mathbf{y}}_{\mathrm{b}\mathrm{u}}^{*}=\mathbf{y}-\mathbf{X}\widehat{\mathbf{b}}-{\sum }_{\begin{array}{l}\mathrm{j}=1\\ \mathrm{j}\ne \mathrm{i}\end{array}}^{\mathrm{f}}{\mathbf{Z}}_{\mathrm{j}}{\widehat{\mathbf{u}}}_{\mathrm{j}}$ = phenotypic observations adjusted for the fixed effects and all random genetic values except those of ${\widehat{\mathbf{u}}}_{\mathrm{i}}$, ${\mathbf{y}}_{\mathrm{u}}^{*}=\mathbf{y}-{\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathbf{Z}}_{\mathrm{i}}{\widehat{\mathbf{u}}}_{\mathrm{i}}$ = phenotypic observations adjusted for all random genetic values, and ${\left({\mathbf{X}}^{\prime }\mathbf{X}\right)}^{-}$ is a generalized inverse of ${\mathbf{X}}^{\prime }\mathbf{X}$. Eq. 31 shows the MME version of ${\widehat{\mathbf{u}}}_{\mathrm{i}}$ uses the phenotypic values adjusted for the GBLUP of all other effect types in the model. Since the MME version of ${\widehat{\mathbf{u}}}_{\mathrm{i}}$ (Eq. 31) and $\widehat{\mathbf{b}}$ (Eq. 32) are identical to the CE version of ${\widehat{\mathbf{u}}}_{\mathrm{i}}$ (Eq. 18) and $\widehat{\mathbf{b}}$ (Eq. 19), the CE version of ${\widehat{\mathbf{u}}}_{\mathrm{i}}$ (Eq. 18) uses the phenotypic values adjusted for the GBLUP of all other effect types in the model even though the CE version does not do such adjustments explicitly.

Genomic restricted maximum estimation (GREML) of variances and heritabilities

The estimation of variance components uses GREML and a combination of EM-REML and AI-REML algorithms of iterative solutions. EM-REML is slow but converges, whereas AI-REML is fast but fails for zero heritability estimates. In our GVCBLUP, GVCHAP, and EPIHAP computing packages that implement these two algorithms (Wang et al., 2014; Prakapenka et al., 2020; Liang et al., 2021), EM-REML is used automatically when AI-REML fails. The EM-REML iterative algorithm for the multifactorial model of Eqs 16 and 17 is

$\begin{array}{cc}{\mathrm{\sigma }}_{\mathrm{i}}^{\mathrm{2(j+1)}}={\mathrm{\sigma }}_{\mathrm{i}}^{\mathrm{2(j)}}\mathbf{y}{\mathbf{P}}^{\left(\mathrm{j}\right)}\mathbf{Z}{\mathbf{S}}_{\mathrm{i}}{\mathbf{Z}}^{\prime }{\mathbf{P}}^{\left(\mathrm{j}\right)}\mathbf{y}/\mathrm{t}\mathrm{r}\left({\mathbf{P}}^{\left(\mathrm{j}\right)}\mathbf{Z}{\mathbf{S}}_{\mathrm{i}}{\mathbf{Z}}^{\prime }\right),& \mathrm{i}=1,\dots ,\mathrm{f}\end{array}(33)$

${\mathrm{\sigma }}_{\mathrm{e}}^{\mathrm{2(j+1)}}={\mathrm{\sigma }}_{\mathrm{e}}^{\mathrm{2(j)}}\mathbf{y}{\mathbf{P}}^{\left(\mathrm{k}\right)}{\mathbf{P}}^{\left(\mathrm{j}\right)}\mathbf{y}/\mathrm{t}\mathrm{r}\left({\mathbf{P}}^{\left(\mathrm{j}\right)}\right)(34)$

where j = iteration number. The AI-REML iterative algorithm is an extension of the early formulations (Johnson and Thompson, 1995; Lee and van der Werf, 2006) to the multifactorial model of Eqs 16 and 17:

${\mathbf{\theta }}^{\left(\mathrm{j}+1\right)}={\mathbf{\theta }}^{\left(\mathrm{j}\right)}+{\left(\mathbf{A}{\mathbf{I}}^{\left(\mathrm{j}\right)}\right)}^{-1}{\Delta }^{\left(\mathrm{j}\right)}(35)$

where $\mathbf{\theta }={\left({\mathrm{\sigma }}_1^2{,\mathrm{\sigma }}_2^2{,...,\mathrm{\sigma }}_{\mathrm{f}}^2,{\mathrm{\sigma }}_{\mathrm{f}+1}^2\right)}^{\prime }$ = $\left(\mathrm{f}+1\right)\times 1$ column vector of variance–covariance components, ${\mathrm{\sigma }}_{\mathrm{f}+1}^2={\mathrm{\sigma }}_{\mathrm{e}}^2$ = residual variance, $\Delta ={\left({\Delta }_1,{\Delta }_2,...,{\Delta }_{\mathrm{f}},{\Delta }_{\mathrm{f}+1}\right)}^{\prime }$ = $\left(\mathrm{f}+1\right)\times 1$ column vector of the partial derivatives of the log residual likelihood function with respect to each variance component, and j = iteration number. A typical term in $\Delta$ (${\Delta }_{\mathrm{i}}$) and a typical term in AI (${\mathrm{A}\mathrm{I}}_{\mathrm{i}\mathrm{k}}$) are

$\begin{array}{l}{\Delta }_{\mathrm{i}}=-\frac{1}{2}\mathrm{t}\mathrm{r}\left(\mathbf{P}\frac{\partial \mathbf{V}}{\partial {\mathrm{\sigma }}_{\mathrm{i}}^2}\right)+\frac{1}{2}{\mathbf{y}}^{\prime }\mathbf{P}\frac{\partial \mathbf{V}}{\partial {\mathrm{\sigma }}_{\mathrm{i}}^2}\mathbf{P}\mathbf{y}\\ =-\frac{1}{2}\mathrm{t}\mathrm{r}\left(\mathbf{P}\mathbf{Z}{\mathbf{S}}_{\mathrm{i}}{\mathbf{Z}}^{\prime }\right)+\frac{1}{2}{\mathbf{y}}^{\prime }\mathbf{P}\mathbf{Z}{\mathbf{S}}_{\mathrm{i}}{\mathbf{Z}}^{\prime }\mathbf{P}\mathbf{y},\mathrm{i}=1,...,\mathrm{f}+1\end{array}(36)$

$\begin{array}{l}{\mathrm{A}\mathrm{I}}_{\mathrm{i}\mathrm{k}}=\frac{1}{2}{\mathbf{y}}^{\prime }\mathbf{P}\frac{\partial \mathbf{V}}{\partial {\mathrm{\sigma }}_{\mathrm{i}}^2}\mathbf{P}\frac{\partial \mathbf{V}}{\partial {\mathrm{\sigma }}_{\mathrm{k}}^2}\mathbf{P}\mathbf{y}\\ \begin{array}{cc}=\frac{1}{2}{\mathbf{y}}^{\prime }\mathbf{P}{\mathbf{Z}}^{\prime }{\mathbf{S}}_{\mathrm{i}}{\mathbf{Z}}^{\prime }\mathbf{P}\mathbf{Z}{\mathbf{S}}_{\mathrm{k}}{\mathbf{Z}}^{\prime }\mathbf{P}\mathbf{y},& \mathrm{i},\mathrm{k}=1,\dots ,\mathrm{f}+1\end{array}\end{array}(37)$

where ${\mathbf{S}}_{\mathrm{f}+1}={\mathbf{I}}_{\mathrm{N}}$. For the full Model-I or Model-II, some effect types inevitably may have zero variances. In those cases, AI-REML (Equations 35–37) fails and EM-REML (Equations 33 and 34) still converges, although a slow convergence rate can be expected for the full Model-I or Model-II. Once the effect types with zero variances are removed from the model, AI-REML converges, and a fast convergence rate can be expected. The estimate of the genomic heritability for each type of genetic effects $\left({\mathrm{h}}_{\mathrm{i}}^2\right)$ and the total heritability of all types of genetic effects $\left({\mathrm{H}}^2\right)$ are

$\begin{array}{cc}{\mathrm{h}}_{\mathrm{i}}^2={\mathrm{\sigma }}_{\mathrm{i}}^2/{\mathrm{\sigma }}_{\mathrm{y}}^2& \mathrm{i}=1,\dots ,\mathrm{f}\end{array}(38)$

${\mathrm{H}}^2={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{h}}_{\mathrm{i}}^2(39)$

where ${\mathrm{\sigma }}_{\mathrm{y}}^2={\sum }_{\mathrm{i}=1}^{\mathrm{f}}{\mathrm{\sigma }}_{\mathrm{i}}^2{+\mathrm{\sigma }}_{\mathrm{e}}^2=$ phenotypic variance.

The heritability estimates of Eq. 38 can be used for model selection by removing effect types with heritability estimates below a user-determined threshold value from the prediction model. Since different traits may have different genetic architectures, we hypothesize that some traits may involve only a small number of the effect types, and some traits are more complex and involve more effect types; global epistasis may be more important than local high-order epistasis effects of haplotypes for some traits, whereas the reverse may be true for other traits, and some traits may be affected by both global high-order and local high-order epistasis effects. The heritability estimates from Eq. 37 provide an approach to evaluate these hypotheses and identify effect types relevant to the phenotypic variance, whereas the total heritability of Eq. 38 provides an estimate of the total genetic contribution to the phenotypic variance. In addition to the use of heritability estimates, prediction accuracy based on GBLUP can be used for model selection by requiring a threshold accuracy level for the effect type to be included in the prediction model, for example, we identified the A + A×A model to have the same accuracy of predicting the phenotypic values of daughter pregnancy rate as the full Model-I in U.S. Holstein cows (Liang et al., 2022).

Estimation of pairwise epistasis effect and heritability

The heritability of an SNP, haplotype block, or pairwise epistasis effect is the contribution of the genetic effect to the phenotypic variance and is also the contribution to the heritability of the effect type, and is estimated through the GBLUP of the corresponding genetic effects. These heritability estimates can be used to identify genome locations with large contributions to the phenotypic variance. The estimation of pairwise epistasis effects and heritability is the most demanding computation because the pairwise epistasis model matrices must be creased and are no longer avoidable. Estimating the effects and heritabilities for third-order epistasis effects is computationally unfeasible and is not considered. GBLUPs of SNP, haplotype, and pairwise epistasis effects of Model-I (Supplementary Table S1) are calculated as

$\begin{array}{cc}{\widehat{\mathbf{\tau }}}_{\mathrm{i}}={\mathrm{\sigma }}_{\mathrm{i}}^2{\mathbf{T}}_{\mathrm{i}}^{\prime }{\mathbf{Z}}^{\prime }\mathbf{P}\mathbf{y}={\mathbf{T}}_{\mathrm{i}}^{\prime }{\mathbf{S}}_{\mathrm{i}}^{-1}{\widehat{\mathbf{u}}}_{\mathrm{i}}& \end{array}(40)$

where ${\widehat{\mathbf{\tau }}}_{\mathrm{i}}$ is the $\mathrm{m}\times 1$ column vector of SNP additive effects for $\mathrm{i}=1$, or SNP dominance effects for $\mathrm{i}=2$; or $\mathrm{b}\times 1$ column vector of haplotype additive effects for $\mathrm{i}=3$; or ${(}_2^{\mathrm{m}})\times 1$ column vector of A×A epistasis effects for $\mathrm{i}=4$, or $2{(}_2^{\mathrm{m}})\times 1$ column vector of A×D epistasis effects for $\mathrm{i}=5$, or ${(}_2^{\mathrm{m}})\times 1$ column vector of D×D epistasis effects for $\mathrm{i}=6$. For $\mathrm{i}=5$; the order of A×D and D×A effects is determined by the order of the model matrices of those effects, that is, ${\widehat{\mathbf{\tau }}}_5={\left({\widehat{\mathbf{\tau }}}_{\mathrm{\alpha }\mathrm{\delta }}^{\prime },{\widehat{\mathbf{\tau }}}_{\mathrm{\delta }\mathrm{\alpha }}^{\prime }\right)}^{\prime }$ if ${\mathbf{T}}_5=\left({\mathbf{T}}_{\mathrm{\alpha }\mathrm{\delta }},{\mathbf{T}}_{\mathrm{\delta }\mathrm{\alpha }}\right)$, or ${\widehat{\mathbf{\tau }}}_5={\left({\widehat{\mathbf{\tau }}}_{\mathrm{\delta }\mathrm{\alpha }}^{\prime },{\widehat{\mathbf{\tau }}}_{\mathrm{\alpha }\mathrm{\delta }}^{\prime }\right)}^{\prime }$ if ${\mathbf{T}}_5=\left({\mathbf{T}}_{\mathrm{\delta }\mathrm{\alpha }},{\mathbf{T}}_{\mathrm{\alpha }\mathrm{\delta }}\right)$. The heritability of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ effect of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type (${\widehat{\mathrm{h}}}_{\mathrm{i}\mathrm{j}}^2$) is estimated as a faction of the genomic heritability of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type (${\widehat{\mathrm{h}}}_{\mathrm{i}}^2$):

${\widehat{\mathrm{h}}}_{\mathrm{i}\mathrm{j}}^2=\left({\widehat{\mathrm{\tau }}}_{\mathrm{i}\mathrm{j}}^2/{\sum }_{\mathrm{i}=1}^{\mathrm{m}}{\widehat{\mathrm{\tau }}}_{\mathrm{i}\mathrm{j}}^2\right){\widehat{\mathrm{h}}}_{\mathrm{i}}^2=\left({\widehat{\mathrm{\tau }}}_{\mathrm{i}\mathrm{j}}^2/{\widehat{\mathbf{\tau }}}_{\mathrm{i}}^{\prime }{\widehat{\mathbf{\tau }}}_{\mathrm{i}}\right){\widehat{\mathrm{h}}}_{\mathrm{i}}^2={\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{y}}^2(41)$

where ${\widehat{\mathbf{\tau }}}_{\mathrm{i}\mathrm{j}}=$ the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ effect of ${\widehat{\mathbf{\tau }}}_{\mathrm{i}},$ ${\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2=$ estimated variance of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type, ${\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2=$ estimated variance of the ${\mathrm{j}}^{\mathrm{t}\mathrm{h}}$ effect of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type, and ${\widehat{\mathrm{h}}}_{\mathrm{i}}^2$ = genomic heritability of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type defined by Equation (52). For proving Equation 57, ${\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2$ and ${\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2$ can be formulated based on the method of mixed-model equations (MME):

${\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2={\widehat{\mathbf{\tau }}}_{\mathrm{i}}^{\prime }{\widehat{\mathbf{\tau }}}_{\mathrm{i}}/\left[{\mathrm{m}}_{\mathrm{i}}-,{\mathrm{t}\mathrm{r}(\mathbf{C}}^{\mathrm{i}\mathrm{i}}{)\mathrm{\lambda }}_{\mathrm{i}}\right]={\sum }_{\mathrm{j}=1}^{{\mathrm{m}}_{\mathrm{i}}}{\mathrm{\tau }}_{\mathrm{i}\mathrm{j}}^2/\left[{\mathrm{m}}_{\mathrm{i}}-{\mathrm{t}\mathrm{r}(\mathbf{C}}^{\mathrm{i}\mathrm{i}}{)\mathrm{\lambda }}_{\mathrm{i}}\right]={\sum }_{\mathrm{j}=1}^{{\mathrm{m}}_{\mathrm{i}}}{\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2(42)$

${{\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2=\widehat{\mathbf{\tau }}}_{\mathrm{i}\mathrm{j}}^{\mathit{2}}/\left[{\mathrm{m}}_{\mathrm{i}}-{\mathrm{t}\mathrm{r}(\mathbf{C}}^{\mathrm{i}\mathrm{i}}{)\mathrm{\lambda }}_{\mathrm{i}}\right],(43)$

where ${\mathbf{C}}^{\mathrm{i}\mathrm{i}}$ is the submatrix in the inverse or generalized inverse of the coefficient matrix of the MME corresponding to the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type, ${\mathrm{m}}_{\mathrm{i}}$ = number of effects of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type, and ${\lambda }_{\mathrm{i}}={\widehat{\mathrm{\sigma }}}_{\mathrm{e}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2$. Dividing Eq. 43 by ${\widehat{\mathrm{\sigma }}}_{\mathrm{y}}^2$ and multiplying by ${\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2$ yield Eq. 41:

${\widehat{\mathrm{h}}}_{\mathrm{i}\mathrm{j}}^2=\left({\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2\right)\left({\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{y}}^2\right)=\left({\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2\right)\left({\widehat{\mathrm{\sigma }}}_{\mathrm{i}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{y}}^2\right)=\left({\widehat{\mathrm{\tau }}}_{\mathrm{i}\mathrm{j}}^2/{\sum }_{\mathrm{i}=1}^{\mathrm{m}}{\widehat{\mathrm{\tau }}}_{\mathrm{i}\mathrm{j}}^2\right){\widehat{\mathrm{h}}}_{\mathrm{i}}^2=\left({\widehat{\mathrm{\tau }}}_{\mathrm{i}\mathrm{j}}^2/{\widehat{\mathbf{\tau }}}_{\mathrm{i}}^{\prime }{\widehat{\mathbf{\tau }}}_{\mathrm{i}}\right){\widehat{\mathrm{h}}}_{\mathrm{i}}^2=\left({\widehat{\mathrm{\sigma }}}_{\mathrm{i}\mathrm{j}}^2/{\widehat{\mathrm{\sigma }}}_{\mathrm{y}}^2\right).$

It is readily seen that the sum of all heritability estimates of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type is the genomic heritability of the ${\mathrm{i}}^{\mathrm{t}\mathrm{h}}$ effect type: ${\sum }_{\mathrm{i}=1}^{{\mathrm{m}}_{\mathrm{i}}}{\widehat{\mathrm{h}}}_{\mathrm{i}\mathrm{j}}^2={\widehat{\mathrm{h}}}_{\mathrm{i}}^2$. It is of note that Eqs 42 and 43 using MME are only for proving Eq. 41. The MME method is computationally prohibitive for estimating genetic effects and their variances under the multifactorial model, although the MME method yields results identical to the CE method, which is computationally feasible for genomic estimation and prediction under the multifactorial model.

Comparison between exact and approximate genomic epistasis relationship matrices

We evaluated the differences between AGERM and EGERM in genomic heritability estimates and prediction accuracies using a publicly available swine genomics data set that had 3,534 animals from a single PIC nucleus pig line with five anonymous traits and 52,842 genotyped and imputed autosome SNPs after filtering by requiring minor allele frequency (MAF) > 0.001 and proportion of missing SNP genotypes < 0.100 (Cleveland et al., 2012). The EGERM followed the method used by Jiang and Reif (2020), and the AGERM methods are described in Supplementary Text S1. The heritability results showed that EGERM had slightly higher heritability estimates than AGERM except for the A×A heritability of T3, where AGERM had a slightly higher estimate than EGERM (0.280 vs. 0.278, Table 1). From Table 1, effect type with nonzero heritability estimates was included in the prediction model for evaluating the observed prediction accuracy as the correlation between the GBLUP of genotypic values and the phenotypic values in each validation population and then averaged over all 10 validation populations. The results showed that AGERM and EGERM had the same prediction accuracy for this swine sample (Table 2). A disadvantage of EGERM is the computing time for the construction of EGERM, about 9.51 times as much time for pairwise relationship matrices, 8.29 as much time for third-order and 9.44 times as much time for fourth-order as required for AGERM (Table 3). However, computing time is not the deciding factor for choosing between the exact and approximate methods because the multi-node approach that calculates each genomic relationship matrix in pieces and adds those pieces together can reduce the computing time to an acceptable level when multiple threads/cores are available, and the two-step strategy can be used so that each genomic relationship is calculated only once for different traits and validation populations (Prakapenka et al., 2020). Prediction accuracy is the ultimate deciding factor in choosing between different methods. We reported results of comparing AGERM and EGERM using 60,671 SNPs and 22,022 first-lactation Holstein cows with phenotypic observations of daughter pregnancy rates, showing that AGERM and EGERM had the same heritability estimates and prediction accuracy, but EGERM required 21 times as much computing time as that required by AGERM, which required 1.32 times as much time for the genomic additive relationship matrix (Liang et al., 2022). The combined results of the swine and Holstein samples indicated that EGERM and AGERM had similar results and that the computing difficulty of EGERM over AGERM increased rapidly as the sample size increased. Given the computing difficulty of EGERM and the negligible differences between EGERM and AGERM in prediction accuracy, AGERM should be favored for its mathematical simplicity and computing efficiency, at least for samples with 50,000 SNPs or more.

TABLE 1

TABLE 1. Genomic heritability estimates of additive, dominance, and epistasis effects up to the third-order for five traits in a swine population.

TABLE 2

TABLE 2. Observed prediction accuracy of epistasis models relative to the additive model for five traits in a swine population.

TABLE 3

TABLE 3. Computing time (in seconds) for the construction of exact and approximate genomic epistasis relationship matrices for a swine population with 3,534 pigs and 52,843 SNPs using 20 threads of the Mangi supercomputer of the Minnesota Supercomputer Institute at the University of Minnesota.

Numerical demonstration

The methods of genomic epistasis relationship matrices based on the additive and dominance model matrices, GREML, GBLUP and reliability, and estimation of effect heritability are demonstrated using an R program (DEMO.R) and a small artificial sample for the convenience of reading the numerical results (Supplementary Text S2 and R program). Because of the artificial nature and the extremely small sample size, this numerical demonstration does not have any genetic and methodology implications and is for showing calculations of the methods only. This R program is an extension of the R demo program of GVCHAP that integrates SNP and haplotype effects and has a computing pipeline for producing the input haplotype data from the SNP data (Prakapenka et al., 2020).

Conclusion

The multifactorial methods with SNP, haplotype, and epistasis effects up to the third-order provide an approach to investigate the contributions of global low-order and local high-order epistasis effects to the phenotypic variance and the accuracy of genomic prediction. Genomic heritability of each effect type from GREML and prediction accuracy from validation studies using GBLUP can be used jointly to identify effect types contributing to the phenotypic variance and the accuracy of genomic prediction, and the GBLUP for the multifactorial model with selected effect type can be used for genomic evaluation. With many capabilities, including the use of intra- and inter-chromosome separately, the multifactorial methods offer a significant methodology capability to investigate and utilize complex genetic mechanisms for genomic prediction and for understanding the complex genome–phenome relationships.

Data availability statement

Publicly available datasets were analyzed in this study. These data can be found at: https://academic.oup.com/g3journal/article/2/4/429/6026060#supplementary-data.

Author contributions

YD conceived this study and derived the formulations. ZL contributed to formulations of the epistasis genomic relationships, implemented the epistasis methods in EPIHAP, and validated and evaluated the methods. DP contributed to the data processing for methodology evaluation. YD and ZL prepared the manuscript.

Funding

This research was supported by the National Institutes of Health’s National Human Genome Research Institute, grant R01HG012425, as part of the NSF/NIH Enabling Discovery through Genomics (EDGE) Program; grant 2020-67015-31133 from the USDA National Institute of Food and Agriculture; and project MIN-16-124 of the Agricultural Experiment Station at the University of Minnesota. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The use of the Ceres and Atlas computers in this research was supported by USDA-ARS projects 8042-31000-002-00-D and 8042-31000-001-00-D.

Acknowledgments

The supercomputer of Minnesota Supercomputer Institute at the University of Minnesota and the Ceres and Atlas high-performance computing system of USDA-ARS were used for the evaluation and testing of the methods and EPIHAP computing package. The authors thank Paul VanRaden, Steven Schroeder, and Ransom Baldwin for help with the use of the USDA-ARS computing facilities.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.922369/full#supplementary-material

References

Bian, C., Prakapenka, D., Tan, C., Yang, R., Zhu, D., Guo, X., et al. (2021). Haplotype genomic prediction of phenotypic values based on chromosome distance and gene boundaries using low-coverage sequencing in Duroc pigs. Genet. Sel. Evol. 53 (1), 78–19. doi:10.1186/s12711-021-00661-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlborg, Ö., and Haley, C. S. (2004). Epistasis: Too often neglected in complex trait studies? Nat. Rev. Genet. 5 (8), 618–625. doi:10.1038/nrg1407

PubMed Abstract | CrossRef Full Text | Google Scholar

Cleveland, M. A., Hickey, J. M., and Forni, S. (2012). A common dataset for genomic analysis of livestock populations. G3 2 (4), 429–435. doi:10.1534/g3.111.001453

PubMed Abstract | CrossRef Full Text | Google Scholar

Cockerham, C. C. (1954). An extension of the concept of partitioning hereditary variance for analysis of covariances among relatives when epistasis is present. Genetics 39 (6), 859–882. doi:10.1093/genetics/39.6.859

PubMed Abstract | CrossRef Full Text | Google Scholar

Cordell, H. J. (2002). Epistasis: What it means, what it doesn't mean, and statistical methods to detect it in humans. Hum. Mol. Genet. 11 (20), 2463–2468. doi:10.1093/hmg/11.20.2463

PubMed Abstract | CrossRef Full Text | Google Scholar

Da, Y. (2015). Multi-allelic haplotype model based on genetic partition for genomic prediction and variance component estimation using SNP markers. BMC Genet. 16 (1), 144. doi:10.1186/s12863-015-0301-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Da, Y., Tan, C., and Parakapenka, D. (2016). 0336 Joint SNP-haplotype analysis for genomic selection based on the invariance property of GBLUP and GREML to duplicate SNPs. J. Animal Sci. 94 (5), 161–162. doi:10.2527/jam2016-0336

CrossRef Full Text | Google Scholar

Da, Y., Wang, C., Wang, S., and Hu, G. (2014). Mixed model methods for genomic prediction and variance component estimation of additive and dominance effects using SNP markers. PLoS One 9 (1), e87666. doi:10.1371/journal.pone.0087666

PubMed Abstract | CrossRef Full Text | Google Scholar

Hayes, B., and Goddard, M. (2010). Genome-wide association and genomic selection in animal breeding. Genome 53 (11), 876–883. doi:10.1139/G10-076

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, C. (1977). Best linear unbiased prediction of breeding values not in the model for records. J. Dairy Sci. 60 (5), 783–787. doi:10.3168/jds.s0022-0302(77)83935-0

CrossRef Full Text | Google Scholar

Henderson, C. (1985). Best linear unbiased prediction of nonadditive genetic merits in noninbred populations. J. Animal Sci. 60 (1), 111–117. doi:10.2527/jas1985.601111x

CrossRef Full Text | Google Scholar

Jiang, Y., and Reif, J. C. (2020). Efficient algorithms for calculating epistatic genomic relationship matrices. Genetics 216 (3), 651–669. doi:10.1534/genetics.120.303459

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., and Reif, J. C. (2015). Modeling epistasis in genomic selection. Genetics 201 (2), 759–768. doi:10.1534/genetics.115.177907

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, D., and Thompson, R. (1995). Restricted maximum likelihood estimation of variance components for univariate animal models using sparse matrix techniques and average information. J. Dairy Sci. 78 (2), 449–456. doi:10.3168/jds.s0022-0302(95)76654-1

CrossRef Full Text | Google Scholar

Kempthorne, O. (1954). The correlation between relatives in a random mating population. Proc. R. Soc. Lond. B Biol. Sci. 143 (910), 102–113.

PubMed Abstract | Google Scholar

Lee, S. H., and van der Werf, J. H. (2006). An efficient variance component approach implementing an average information REML suitable for combined LD and linkage mapping with a general complex pedigree. Genet. Sel. Evol. 38 (1), 25–43. doi:10.1186/1297-9686-38-1-25

PubMed Abstract | CrossRef Full Text | Google Scholar

Liang, Z., Prakapenka, D., and Da, Y. (2022). Comparison of two methods of genomic epistasis relationship matrices using daughter pregnancy rate in U.S. Holstein cattle. Abstract 2466V, page 409 of ADSA2022 Abstracts. Available at: https://www.adsa.org/Portals/0/SiteContent/Docs/Meetings/2022ADSA/Abstracts_BOOK_2022.pdf?v=20220613 (Accessed September 27, 2022).

Google Scholar

Liang, Z., Prakapenka, D., and Da, Y. (2021). Epihap: A computing tool for genomic estimation and prediction using global epistasis effects and haplotype effects. Abstract P167, page 223 of ADSA2021 Abstracts, ADSA 2021 Virtual Annual Meeting. Available at: https://www.adsa.org/Portals/0/SiteContent/Docs/Meetings/2021ADSA/ADSA2021_Abstracts.pdf (Accessed September 27, 2022).

Google Scholar

Liang, Z., Tan, C., Prakapenka, D., Ma, L., and Da, Y. (2020). Haplotype analysis of genomic prediction using structural and functional genomic information for seven human phenotypes. Front. Genet. 11 (1461), 588907. doi:10.3389/fgene.2020.588907

PubMed Abstract | CrossRef Full Text | Google Scholar

Mackay, T. F. (2014). Epistasis and quantitative traits: Using model organisms to study gene–gene interactions. Nat. Rev. Genet. 15 (1), 22–33. doi:10.1038/nrg3627

PubMed Abstract | CrossRef Full Text | Google Scholar

Martini, J. W., Toledo, F. H., and Crossa, J. (2020). On the approximation of interaction effect models by Hadamard powers of the additive genomic relationship. Theor. Popul. Biol. 132, 16–23. doi:10.1016/j.tpb.2020.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Martini, J. W., Wimmer, V., Erbe, M., and Simianer, H. (2016). Epistasis and covariance: How gene interaction translates into genomic relationship. Theor. Appl. Genet. 129 (5), 963–976. doi:10.1007/s00122-016-2675-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz, P. R., Resende, M. F., Gezan, S. A., Resende, M. D. V., de los Campos, G., Kirst, M., et al. (2014). Unraveling additive from nonadditive effects using genomic relationship matrices. Genetics 198 (4), 1759–1768. doi:10.1534/genetics.114.171322

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, P. C. (2008). Epistasis—The essential role of gene interactions in the structure and evolution of genetic systems. Nat. Rev. Genet. 9 (11), 855–867. doi:10.1038/nrg2452

PubMed Abstract | CrossRef Full Text | Google Scholar

Prakapenka, D., Liang, Z., Jiang, J., Ma, L., and Da, Y. (2021). A Large-scale genome-wide association study of epistasis effects of production traits and daughter pregnancy rate in US Holstein cattle. Genes 12 (7), 1089. doi:10.3390/genes12071089

PubMed Abstract | CrossRef Full Text | Google Scholar

Prakapenka, D., Wang, C., Liang, Z., Bian, C., Tan, C., and Da, Y. (2020). Gvchap: A computing pipeline for genomic prediction and variance component estimation using haplotypes and SNP markers. Front. Genet. 11, 282. doi:10.3389/fgene.2020.00282

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. D., and Van Steen, K. (2018). The search for gene-gene interactions in genome-wide association studies: Challenges in abundance of methods, practical considerations, and biological interpretation. Ann. Transl. Med. 6 (8), 157. doi:10.21037/atm.2018.04.05

PubMed Abstract | CrossRef Full Text | Google Scholar

Segre, D., DeLuna, A., Church, G. M., and Kishony, R. (2005). Modular epistasis in yeast metabolism. Nat. Genet. 37 (1), 77–83. doi:10.1038/ng1489

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, G., Christensen, O. F., Ostersen, T., Henryon, M., and Lund, M. S. (2012). Estimating additive and non-additive genetic variances and predicting genetic merits using genome-wide dense single nucleotide polymorphism markers. PloS One 7 (9), e45293. doi:10.1371/journal.pone.0045293

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, C., Wu, Z., Ren, J., Huang, Z., Liu, D., He, X., et al. (2017). Genome-wide association study and accuracy of genomic prediction for teat number in Duroc pigs using genotyping-by-sequencing. Genet. Sel. Evol. 49 (1), 35. doi:10.1186/s12711-017-0311-8

PubMed Abstract | CrossRef Full Text | Google Scholar

VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. J. Dairy Sci. 91 (11), 4414–4423. doi:10.3168/jds.2007-0980

PubMed Abstract | CrossRef Full Text | Google Scholar

Vitezica, Z. G., Legarra, A., Toro, M. A., and Varona, L. (2017). Orthogonal estimates of variances for additive, dominance, and epistatic effects in populations. Genetics 206 (3), 1297–1307. doi:10.1534/genetics.116.199406

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., and Da, Y. (2014). Quantitative genetics model as the unifying model for defining genomic relationship and inbreeding coefficient. PLoS One 9 (12), e114484. doi:10.1371/journal.pone.0114484

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, C., Prakapenka, D., Wang, S., Pulugurta, S., Runesha, H. B., and Da, Y. (2014). Gvcblup: A computer package for genomic prediction and variance component estimation of additive and dominance effects. BMC Bioinforma. 15 (1), 270. doi:10.1186/1471-2105-15-270

PubMed Abstract | CrossRef Full Text | Google Scholar