Reduced Animal Models Fitting Only Equations for Phenotyped Animals

Reduced models are equivalent models to the full model that enable reduction in the computational demand for solving the problem, here, mixed model equations for estimating breeding values of selection candidates. Since phenotyped animals provide data to the model, the aim of this study was to reduce animal models to those equations corresponding to phenotyped animals. Non-phenotyped ancestral animals have normally been included in analyses as they facilitate formation of the inverse numerator relationship matrix. However, a reduced model can exclude those animals and obtain identical solutions for the breeding values of the animals of interest. Solutions corresponding to non-phenotyped animals can be back-solved from the solutions of phenotyped animals and specific blocks of the inverted relationship matrix. This idea was extended to other forms of animal model and the results from each reduced model (and back-solving) were identical to the results from the corresponding full model. Previous studies have been mainly focused on reduced animal models that absorb equations corresponding to non-parents and solve equations only for parents of phenotyped animals. These two types of reduced animal model can be combined to formulate only equations corresponding to phenotyped parents of phenotyped progeny.


INTRODUCTION
Computational limitations have been a major challenge facing animal breeders, especially for solving large mixed model equations, such as those including multiple traits and millions of animals. Computational speed, power and technology have advanced, making most of the computations that were impossible in the past, feasible. However, the amount of data is growing rapidly, and the complexity of evaluation models are increasing. Henderson (1974) introduced the concepts of equivalent models and reduced models to animal breeding. Any linear model can be written in various forms of equivalent models yielding the same first and second moments of the data (Henderson, 1985). A reduced model is an equivalent model to the full model, allowing computational simplifications by reduction in the number of equations to be solved.
Limited by available memory to store the entire set of equations, Henderson (1974) proposed an equivalent reduced model in which the equations for fixed herd effects were absorbed into the equations for sire effects. The full model explicitly fitted herd and sire effects. Generally, absorption is a model rank reduction (reduced number of equations to be solved), by which some effects get absorbed into other effects with no influence on the solutions of the remaining effects in the model.
Reduced models that absorbed multiple-trait equations for random effects were first introduced to animal models by Quaas and Pollak (1980), where they showed that the best linear unbiased predictions (BLUP) of parents (with phenotyped progeny) can be obtained from a reduced model that only includes equations for breeding values for parents. This was achieved by absorbing equations of non-parent individuals (and parents with no phenotyped progeny). The BLUP of breeding values for nonparent individuals can then be back-solved given the BLUP of breeding values for their parents. The concept of reduced models was extended to animal models incorporating quantitative trait loci, introduced by Fernando and Grossman (1989), and was also based on absorbing equations corresponding to non-parent individuals (Cantet and Smith, 1991;Hoeschele, 1993;Arendonk et al., 1994;Iwaisaki, 1996, 1997). Henderson (1986) proposed the estimation of variance components in reduced (based on the methodology of Quaas and Pollak, 1980) rather than full animal models for single traits and single records, based on minimum variance quadratic unbiased estimation (MIVQUE), restricted maximum likelihood (REML) by iterated MIVQUE, and REML by the expectationmaximization (EM, Dempster et al., 1977) algorithm. Providing numerical examples, he showed that the reduced model can be used for estimation of variances, converging to the same solutions obtained by the full model.
This study revisited the reduced animal model from a different perspective by only including equations corresponding to phenotyped animals in the model. In a full model, non-phenotyped parents do not contribute any phenotypic information to the model, but make it easier to account for similarity between animals with phenotypes when they are included in the construction of the inverse numerator relationship matrix. This form of reduced animal model that absorbs the effects of non-phenotyped breeding values is first explained in the context of a simple univariate animal model, and then extended to several other forms of animal model. Numerical examples and R scripts are provided to facilitate understanding.

METHODS
Animal models can be expressed and fitted in various forms and combinations. In this section, some equivalent forms of animal model, their size reduction (absorption) to phenotyped animals (denoted as p), and approach to back-solve solutions corresponding to non-phenotyped animals (denoted as n) are presented.

Basic Animal Model
In its simplest commonly-used form (y = Xb + Zu + e), an animal model contains the vectors of phenotypes (y), fixed effects (b), random direct additive genetic effects [u ∼ N(0, Aσ 2 u ), representing breeding values] and random residuals (e ∼ N(0, Iσ 2 e )), with σ 2 u and σ 2 e being the corresponding known variances, and A being the numerator relationship matrix. Matrices X and Z relate phenotypes to fixed effects and animals, respectively. In matrix notation solutions to this single-trait model are obtained by solving the equation system: whereb is a vector of solutions for fixed effects (b) andû is the vector of best linear unbiased predictions for u (known as estimated breeding values) and λ = σ 2 e /σ 2 u . Permuting rows and columns of the sparse matrix A −1 so that non-phenotyped animals appear before phenotyped animals, where 0 is a matrix with elements equal to zero. Here, 0 has number of rows equal to the number of phenotypes and number of columns equal to the number of non-phenotyped animals. Then, Equation (1) can be written as: Permuting y so that animals are in the same order as A pp , Z p = I. Then, If we absorb the equations corresponding to breeding values for non-phenotyped animals by modifying the equations corresponding to breeding values for phenotyped animals, Equation (3) is transformed to Equation (4), which contains breeding values only for phenotyped animals: Derivation of is provided later in this subsection. Considering equations corresponding to breeding values for non-phenotyped animals in Equation (3), we have the identity: Thus,û n can be back-solved from A nn , A np , andû p . Comparing Equations (1) and (4), y = Xb + Zu + e is reduced to y = Xb + u p + e. Thus, phenotypes depend only on the genetic merit of phenotyped animals and are independent from the genetic merit of their non-phenotyped relatives. The reduced model has a different structure of additive genetic variance (u p ∼ N(0, −1 σ 2 u ) vs. u ∼ N(0, Aσ 2 u )). Consider the typically dense matrix A = A nn A np A pn A pp . The following mixed model equations are equivalent to Equation (1), Frontiers in Genetics | www.frontiersin.org but can only be used to directly predictb andû p , and it involves no changes in model assumptions (Henderson, 1986): The above representation of the model is useful to prove the derivation of , but not computationally efficient, because A pp is generally not sparse and it cannot be inverted using the indirect inversion methods of Henderson (1976) or Quaas (1976). Equating row 3 of Equation (3) (Xb+λA pnû n +(I+λA pp )û p = y) and row 2 of Equation (7) Solving the obtained equation with Equation (6): Thus, forming and solving the mixed model equations using will produce the same results as forming and solving the mixed model equations using A −1 pp , but without the need to form A, invert or store A pp .
We carried out a test to compare the computational time between the full and the reduced animal model. Pedigree and phenotypes of 10,000 animals were simulated using R package pedSimulate (Nilforooshan, 2021). The pedigree was initiated with 100 founder animals (generation 0) of both sexes (50 each) and continued to eight generations. Except 2,000 randomly selected phenotypes from generations 6 to 8, other phenotypes were set to missing. Phenotypes were randomly assigned to five herds and the herd effect was added to phenotypes. The R statistical software (R Core Team, 2019) was used to analyse the data.

Multiple Trait Animal Model
The mixed model equations for a multiple trait animal model in its simplest form with no additional random effects other than for breeding values can be represented as: where G is the t × t matrix of genetic covariances between the traits, R is the matrix of residual covariances between all the residuals on all the phenotyped animals, and ⊗ is the Kronecker product (Horn and Johnson, 1991). Permuting rows and columns of A −1 so that animals that are non-phenotyped (for all the traits) appear before animals that are phenotyped for at least one of the traits, Z = 0 Z p , Z p corresponds to animals phenotyped for at least one of the traits, and: In the case of a 2-trait model, X = X 1 0 0 X 2 and Z p = Z p 1 0 0 Z p 2 . Z p 1 and Z p 2 contain columns for the same set of animals (those phenotyped for either or both of the two traits). Following the same procedure as for the single-trait animal model, the model is reduced to fitting phenotyped animals for at least one of the traits in the model: Considering Equation (9), solutions for non-phenotyped animals (û n ) can be back-solved as:

Repeatability Animal Model
Repeatability animal models are employed for the analysis of repeated phenotypes of the same trait on an animal. In the simplest repeatability animal models, a genetic correlation of 1 and equal residual correlations are considered between repeated records from an animal, and the same variances are applied to all the phenotypes (Mrode, 2005). Repeatability animal models are typically fitted including an additional random term for permanent environmental effects (y = Xb + Zu + Wp + e). In matrix notation: wherep is the vector of solutions for random permanent environmental effects, W is the incidence matrix relating phenotypes to permanent environmental effects), α = σ 2 e /σ 2 p , and σ 2 p is the permanent environment variance. Considering W being already limited to phenotyped animals, W = Z p , and: (13) Equation (13) can be reduced to phenotyped animals (Equation 14) by absorbing the equations corresponding to breeding values for non-phenotyped animals into the equations corresponding to breeding values for phenotyped animals: Back-solving solutions for non-phenotyped animals is done according to Equation (6). If permanent environmental effect solutions are not of interest, Equation (14) can be further reduced by absorbing the equations for random permanent environmental effects. Any random effect not correlated to other random effects remaining in the model can be confounded with the residual term (Quaas and Pollak, 1980).

Maternal Trait Animal Models
For some traits, dams not only influence the phenotypic expression of their progeny via directly passed genes, but also through maternal ability that influences the environment for their progeny, such as the milk volume they provide to the offspring (Mrode, 2005). Maternal ability is partly genetic and partly environmental. Accordingly, there can be animal models with maternal genetic effects, maternal permanent environment effects, or both. Animal models with maternal permanent environment effects have similar structure to repeatability animal models with the differences that there might be no repeated records per animal, W is the incidence matrix relating phenotypes to the dams of phenotyped animals, and a different variance than σ 2 p is involved (σ 2 mpe , maternal permanent environment variance). Maternal permanent environmental effect is already limited to phenotyped animals as the incidence matrix used for this effect has rows corresponding to phenotypes and columns corresponding to dams of phenotyped animals.
Animal models with maternal genetic effects (y = Xb + Zu + Sm + e) uncorrelated to direct additive genetic effects have the variance structure: where σ 2 m is the maternal genetic variance. In matrix notation: Equation (15) can be re-written as: where λ 1 and λ 2 are equal to σ 2 e /σ 2 u and σ 2 e /σ 2 m , respectively. Considering equations corresponding toû n in Equation (16), Equation (6) is derived. Equation (16) can then be reduced to: Equations corresponding tom cannot be reduced to equations corresponding tom p . The reason is the non-zero S n matrix. Animal models with maternal genetic effects correlated to direct additive genetic effects have the variance structure: In matrix notation: For an efficient model reduction, it is important that all A −1 occurrences in Equation (18) are converted to the same matrix. In order to do that, the definition of n and p are changed to p being phenotyped animals or non-phenotyped dams, and n being other animals. That enables both Z n and S n being zero matrices. Then, permuting rows and columns of A −1 so that animals n appear before animals p, Equation (18) can be written as: Reducing the above equation to phenotyped animals, equations corresponding toû n are absorbed in equations corresponding tô u p , and equations corresponding tom n are absorbed in equations corresponding tom p , resulting to: u n andm n are back-solved givenû p ,m p , A nn , A np , and G −1 . For the proof, consider rows 2 and 4 in Equation (19):

Animal Models With Genetic Groups
The idea of sire models with genetic groups (Quaas and Pollak, 1981) was extended to animal models by Westell and Van Vleck (1987), where grouping stems from the fact that unknown parents belonging to different genetic backgrounds and periods of time have different genetic means. An animal model with genetic groups and no additional random effect other than breeding values is written as y = Xb + Zu + ZQg + e. In matrix notation: whereĝ is the vector of solutions for genetic groups and Q is the matrix relating genetic groups to animals. Permuting rows and columns of A −1 so that non-phenotyped animals appear before phenotyped animals and phenotyped animals appear in the same order in y and A pp , Z = 0 Z p , Q ′ = Q ′ n Q ′ p , ZQ = Z p Q p , and Z p = I. Thus, Equation (22) can be written as: Absorbing the equations corresponding toû n in the equations corresponding toû p , the above model is reduced to: Back-solvingû n is done using Equation (6) and the breeding value of animals are calculated asû + Qĝ. Forming Q to calculate breeding values outside of the animal model is an additional computational step. To avoid that step, Quaas and Pollak (1981) transformed the mixed model equations so that the solution for the breeding value includes the genetic group effect: Let's consider: where g corresponds to genetic groups in a larger numerator relationship matrix (i.e., is larger than A). Then, Equation (25) can be written as: Considering Z pg = Z p 0 , in which 0 is a matrix of zeros with columns corresponding to genetic groups: Frontiers in Genetics | www.frontiersin.org Applying the same technique used to reduce Equation (3) to Equation (4), Equation (26) is reduced to Equation (27): where extending Equation (5) to include genetic groups: and extending Equation (6) to include genetic groups:

Animal Models With a Genomic Relationship Matrix (GBLUP)
In the previous sections, animal models with pedigree-based relationship matrices were discussed. With the availability of genomic data, various forms of genomic relationship matrix (e.g., VanRaden, 2008;Hickey et al., 2013) are used in animal models.
Replacing A with the genomic relationship matrix G in Equation (1) results in: However, it should be noted that unlike Equation (1), Equation (30) is limited to genotyped animals. Due to efficient methods for deriving A −1 without calculating and inverting A (Henderson, 1976;Quaas, 1976), Equation (4) is preferred over Equation (7). However, those methods are not applicable to G, and this form of the equations requires G to be calculated and directly inverted. As a result, the reduced method of choice is equivalent to Equation (7): Therefore, only G pp needs to be formed and inverted. If the calculation of G pp requires allele frequencies, estimation of allele frequencies should not be limited to phenotyped animals. Then, u n can be derived indirectly by, first back-solving marker effects (â) fromû p (Strandén and Garrick, 2009), and then summing up marker effect solutions over all markers for non-phenotyped animals (û n = M nâ (Meuwissen et al., 2001), where M n contains rows of M for non-phenotyped animals, and M is the genotype matrix with coefficients −1, 0, and 1).

Animal Models With an Augmented
Pedigree-Genomic Relationship Matrix (ssGBLUP) Aguilar et al. (2010) and Christensen and Lund (2010) introduced an animal model for simultaneous genetic evaluation of genotyped and non-genotyped animal, using a pedigree-genomic augmented inverse matrix (H −1 ) to be used instead of A −1 in the mixed model equations, where H 11 corresponds to the block for non-genotyped animals, and H 22 corresponds to the block for genotyped animals. Though, it is possible to partition H −1 to H nn H np H pn H pp and reduce the model to phenotyped animals (similar to the procedures already explained), perhaps a more practical approach would be to split animals into three groups of "non-phenotyped and non-genotyped" animals (4), "phenotyped and non-genotyped" animals (3), and genotyped animals (2) Grouping phenotyped or genotyped animals together (5): Then, an animal model like the following: can be reduced to: where = H 55 − A 54 (A 44 ) −1 A 45 and Z 5 contains columns of Z corresponding to group 5 of animals. After obtainingû 5 ,û 4 can be back-solved by solving

RESULTS
Results from numerical examples are provided. The numerical examples were chosen from small datasets used in Mrode (2005) for different animal models. R scripts for producing the data and the analyses are available in the data repository. Files for any given filename in this study are also available in the data repository. Back-solving were done by direct solving of Equations (6), (11), (21), and (29). Where solving these equations for large pedigree is computationally challenging, preconditioned conjugate gradient method (Strandén and Lidauer, 1999) can be adopted.

Basic Animal Model
Data from Example 3.1 of Mrode (2005) was used (Supplementary Table 1). The data can be produced using e31data.R in the data repository. Equation (1) is coded in function "BLUP" in e31.R. The following code produces the data and solutions.
# Create data source("e31data.R") # Call the solver function source("e31.R") # Solve the full model sol1 = BLUP(X, y, Z, lambda, Ai) Applying Equation (4)  Following Equation (6),û n are back-solved using the following code: u_n = solve(mats$Ann, mats$Anp % * % -sol2[-2:-1]) The following command confirms that solutions from the full (Equation 1) and the reduced models (Equation 4) are identical. We run a test to compare the computational time for the full and the reduced animal model for a pedigree of 10,000 animals, of which 2,000 animals had been phenotyped. The runtime were 10 and 8.8 s for the full and the reduced animal model, respectively. Additionally, it took 2 s for the reduced animal model to calculate . Matrix was considerably denser than the A −1 matrix.

Multiple Trait Animal Model
Data from Example 5.3 of Mrode (2005) was used (Supplementary Table 2). The data can be produced using e53data.R in the data repository. Equation (8) is coded in function "BLUP" in e53.R. The following code produces the data and solutions.

Repeatability Animal Model
Data from Example 4.1 of Mrode (2005) was used (Supplementary Table 3). The data can be produced using e41data.R in the data repository. Equation (12) is coded in function "BLUP" in e41.R. The following code produces the data and solutions.

Maternal Trait Animal Models
Data from Example 6.1 of Mrode (2005) was used (Supplementary Table 4). This example includes maternal genetic and maternal permanent environment effects. The data can be produced using e61data.R in the data repository. The equation system is coded in function "BLUP" in e61.R. The following code produces the data and solutions. The following command confirms that solutions from the full and the reduced models are identical.

Animal Models With Genetic Groups
Data from Example 3.4 of Mrode (2005) was used (Supplementary Table 5). The data can be produced using e34data.R in the data repository. Equation (25) is coded in function "BLUP" in e34.R. The following code produces the data and solutions.
# Create data source("e34data.R") # Call the solver function source("e34.R") # Solve the full model sol9 = BLUP(X, y, Z, lambda, Ai, Q) Applying Equation (27) instead of Equation (25) requires data transformation. Matrices , A nn , and A np are obtained using function "getPhi" coded in getPhi.R. Solutions of the reduced model are obtained using function "BLUP" coded in e31.R, and the following code.

Animal Models With a Genomic Relationship Matrix (GBLUP)
Data from Example 3.1 of Mrode (2005) was used (Supplementary Table 1). The data can be produced using e31data.R in the data repository. Rather than using the A −1 matrix produced by e31data.R (data repository), a G −1 matrix was used, created by makeGi.R (data repository). The following code creates the data and solutions. The solutions are obtained using the function "BLUP" in e31.R.

DISCUSSION
Genetic evaluation centers around the world deal with tight deadlines for delivering breeding value predictions of animals to the industry. The amount of data to be dealt with usually contains millions of animals and that can make restraints around computational time and resources. Often, there are many historical or ancestral animals as well as inactive (culled) breeding animals, for which breeding values are not routinely needed. Computational demand can be reduced by applying reduced animal models that are equivalent to the full animal models.
Several forms of animal model were reduced to comprise equations for breeding values only for phenotyped animals, and extension to other forms of animal model can be done following the same logic. The efficiency of these models depends on several factors, such as the proportion of phenotyped animals in the population being evaluated, sparseness of , and the number of effects in common between the full and the reduced model. There is also a trade-off between the number of phenotyped and non-phenotyped animals. The smaller the proportion of phenotyped animals the smaller the size of MME, but the greater the computational cost of calculating as a result of a larger A nn . The sparsity of , which itself depends on the sparsity of A nn [the sparser the A nn , the sparser the (A nn ) −1 ] can heavily reduce the efficiency of the reduced animal model. Depending on the size of , it might have more non-zero elements than the A −1 matrix. This might erode any computational benefit from the reduced model compared to the full model. We suggest two solutions for this problem. One is discarding non-phenotyped animals that have little or no relatedness with phenotyped progeny (i.e., a few generations distant from any phenotyped descendant), and the other is re-ordering A nn to obtain two diagonal blocks, one very sparse (A ss ) and the other (A dd ) relatively dense Then, A dd , A sd , and A ds are appended to A pp , A np , and A pn . This also requires appending rows of 0 to X and columns of 0 to Z (corresponding to the appended nonphenotyped animals). A pedigree-only-based back-solving (such as for BLUP and for ssGBLUP dividing animals to "phenotyped or genotyped" and "non-phenotyped and non-genotyped") can be done efficiently for millions of animals.
Solutions for non-phenotyped animals can be back-solved after the solutions for phenotyped animals are obtained. Generally, back-solving is more efficient and faster than solving the mixed model equations, and that is due to the known structure of A −1 . Usually, breeding values of only a small proportion of non-phenotyped animals are of interest (e.g., active non-phenotyped breeding animals). Consequently, the back-solving procedure can be omitted by explicitly including those animals together with phenotyped animals in the reduced animal model. The same reduction technique is applicable to BLUP using the genomic relationship matrix (GBLUP). That would reduce the computational demand to form and invert G by limiting it to phenotyped animals only, and the computational demand to solve a smaller set of mixed model equations. Note that the computational cost of inverting G is cubic to its dimension unless it is not full rank because it was formed with fewer markers than the number of genotyped animals, where an equivalent model (Fernando et al., 2016) can be used or it uses a reduced rank approximation, such as the algorithms introduced by Misztal (2016) and Mäntysaari et al. (2017). Though, back-solving is possible, it is discouraged as it involves the availability of the blocks of G −1 among non-phenotyped animals and between non-phenotyped and phenotyped animals. Because most of the computational complexity of GBLUP is due to forming and inverting G, it is recommended to form and invert G for phenotyped animals together with non-phenotyped animals of interest. However, it is assumed that either the allele frequencies in the whole population are known or allele frequencies in the individuals present in the reduced model are not different from the allele frequencies in the whole population. Also, as it was mentioned earlier, another way of back-solving for nonphenotyped animals is through back-solving of marker effects. There are other forms of reduced GBLUP proposed by other studies (Cantet and Smith, 1991;Hoeschele, 1993;Arendonk et al., 1994;Iwaisaki, 1996, 1997). Comparing the efficiency of these methods and possibly combining them with the reduced GBLUP presented in this study are yet to be investigated. The computational efficiency of the reduced GBLUP presented in this study heavily depends on the proportion of animals not included in the reduced model (genotyped and non-phenotyped animals whose evaluation is not of interest).
There are two possibilities for ssGBLUP; dividing animals to phenotyped and non-phenotyped animals or dividing animals to "phenotyped or genotyped" and "non-phenotyped and nongenotyped" groups. The latter is preferred as it makes backsolving independent from H 22 (block of H −1 corresponding to genotyped animals) and dependent to blocks of A −1 between the two groups of animals and among "non-phenotyped and nongenotyped" animals. Blocks of A −1 are much easier to calculate and keep in the memory, because of its sparsity.
Variance components can be estimated using reduced animal models (e.g., Henderson, 1986;Besbes et al., 1992;White et al., 2006;Guy et al., 2009). Restricted maximum likelihood (REML, Patterson and Thompson, 1971) has been a popular method for variance components estimation. Generalized inverses of coefficient matrices of mixed model equations are required under both the full and the reduced animal models (Henderson, 1986). Reduced animal models have considerably smaller order coefficient matrix and consequently less computational complexity of variance components estimation. Henderson (1986) estimated variances for the reduced model based on the method of Quaas and Pollak (1980). Parents with phenotyped progeny are evaluated in this reduced model. To account for Mendelian sampling and for the genetic contribution of unknown parents, there is an additional component (compared to the full model) to the variance of residuals (Henderson, 1986).
The introduced reduction method in this study evaluates all the phenotyped animals and no extra assumption or component is required to hold the assumption on var(y). The matrices used in the reduced model do not make any deviation from the assumptions on variance components. Also, creating the Z matrix for the reduced model in this study is trivial (if the vector of phenotypes is ordered by the order of animals in the pedigree, Z is an identity matrix with the order of the number of phenotypes; for a repeatability model, it is still an incidence matrix with coefficients 0 and 1, relating phenotypes to animals). Creating the Z matrix for the reduced model of Quaas and Pollak (1980) requires pedigree information.
The presented reduced animal models can further be reduced. This can be achieved by a "2-fold" model reduction, which involves reducing the model to equations for phenotyped animals followed by reducing the model to parents of phenotyped animals (Quaas and Pollak, 1980). Such 2-fold reduced animal model contains equations corresponding to phenotyped parents of phenotyped animals. A 2-fold reduced animal model also requires a 2-fold back-solving procedure, in which the solutions of non-phenotyped parents are back-solved from the solutions of phenotyped parents, and then the solutions of non-parents (and parents of non-phenotyped progeny) are back-solved from the solutions of parents. This 2-fold back-solving can be omitted or reduced to a 1-fold back-solving by putting together a group of animals of interest to phenotyped parents with phenotyped progeny. Further research is required on this subject.
Finally, depending on the size of the data and the number of traits and effects in the model, the equation system of the reduced animal model may become small enough that direct solving of the mixed model equations might become faster and more efficient than iterative procedures for solving mixed model equations (Schaeffer and Kennedy, 1986;Berger et al., 1989;Strandén and Lidauer, 1999).

DATA AVAILABILITY STATEMENT
The R scripts to produce the data and perform the analyses for this study can be found in the figshare repository [https://doi.org/ 10.6084/m9.figshare.13369607].

AUTHOR CONTRIBUTIONS
DG initiated the idea of the study. MN extended the idea to other forms of animal model, drafted the manuscript, coded the examples, and validated the results. Both DG and MN took active part in the preparation and revision of the manuscript.

FUNDING
This study received financial support from the NZ Ministry of Primary Industries, SFF Futures Programme: Resilient Dairy-Innovative breeding for a sustainable dairy future (Grant number: PGP06-17006).