Maize populations 1 and 2

In each one of the two F₂ maize (Zea mays) Populations, three variables were recorded: grain yield (RG, Mg ha^‒1), ear height (AM, cm), and plant height (AP, cm). The maize Population 1 had 199 MM and 247 genotypes, whereas in the maize Population 2 the number of MM was 259 and of genotypes 248. The estimated genetic correlations between RG and AM, RG and AP, and AM and AP in the maize Population 1 were, respectively, 0.53, 0.52 and 0.98, whereas in the maize Population 2, those correlations were: 0.58, 0.76 and 0.71.

Population 3 (wheat population)

The wheat (Triticum aestivum L.) double haploid Population included 1279 MM and 599 genotypes. In it, the grain yield (RG, Mg ha^‒1) was recorded in three environments (RG1, RG2, RG3). To predict the breeding value of the candidates for selection, RG1, RG2 and RG3, each of them was considered a particular characteristic because the genotypes were evaluated in different environments. The estimated genetic correlations between RG1 and RG2, RG1 and RG3, and RG2 and RG3, were ‒0.03, ‒0.21 and 0.73, respectively.

The proposed univariate linear model

Let γq=Xuq be a g×1 vector (g = number of genotypes in the population) of genomic breeding values associated to the characteristic q(q=1,2,...,t; t=number of variables) of the candidates for selection. Assume that 𝛾 _q has multivariate normal distribution (NMV) with mean 0 and variance Gσγq2, i.e., 𝛾 _q ~ NMV 0,Gσγq2, where σγq2 is the additive genomic variance of 𝛾 _q and G = XX′ / k is the g×g additive genomic relationships matrix between genotypes; X is a g× m matrix (m= number of MM in the population) of coded MM values (2‒2p, 1‒2p and ‒2p for the genotypes AA, Aa, and aa, respectively) associated to the additive effects of the quantitative traits loci (QTL); p is the frequency of the allele A and 1‒p is the frequency of allele a in the MM j (j=1,2,..., m); u _q is a m ×1 vector of additive effects of the QTL associated to the m MM that affect trait q; k=∑j=1m2pj1-pj (^{Habier et al., 2007}) in a population F₂ and k=∑j=1m4pj1-pj in a double haploid population. In addition, let a _q ~ NMV 0,Aσaq2 be a g×1 vector of additive genetic merits unexplained by the MM associated to the trait q, where A is the numerical relationships matrix and σaq2 is the additive genetic variance of a _q . The combined linear model for the trait qyq* can be denoted as yq*=1μq+Zaq+Zγq+eq, or equivalently as:

(1)

where y _q = yq* −1μ _q ~ NMV 0,Vq is a vector g×1 of the q trait observations centered with regard to the q trait mean, 𝛍 _q ; 1 is a g×1 vector of ones; Vq=Aσaq2+Gσγq2+2Covaq,γq´+Igσeq2, and Covaq,γq´=Gσγq2 (i.e., the covariance between a _q and 𝛾 _q is equal to the variance of 𝛾 _q ); 𝛾 _q , G y σγq2 were defined before; Z is a matrix of incidence (generally an identity matrix g×g) y e _q ~ NMV 0,Igσeq2 is a g×1 vector of residuals; I _g is an identity matrix g×g and σeq2 is the variance of residuals. The model of Equation 1 will be called univariate combined linear model.

Posterior joint distribution of a _q and 𝛾 _q

The posterior joint distribution of a _q and 𝛾 _q can be written as:

(2)

where the symbol “∝” indicates that P(a _q , 𝛾 _{q /} y _q ) can be written as the product of the likelihood function, of y _q , P(y _q / a _q /𝛾 _q ) ∝ exp-12yq-Zaq-Zγq⁡´ R-1yq-Zaq-Zγq, the conditional a priori distribution of a _q given 𝛾 _q , P(a _q / 𝛾 _q ) ∝ exp-12aq-γq´T-1aq-γq and the a priori distribution of 𝛾 _q , P(𝛾 _q ) ∝ exp-12γq´Φ-1γq, R=Igσeq2, T=Aσaq2-Φ y Φ=Gσγq2. According to the properties of the NMV distribution (^{Sorensen and Gianola, 2002}), 𝛾 _q y T are the expectation and the variance of a _q / 𝛾 _q , respectively. Thus, Equation 2 is equal to:

(3)

The right side of Equation 3 is the normal distribution kernel with mean Dd and variance D, where θq´=aq´γq´, D-1=D11-1D12-1D21-1D22-1^-1, D11-1=R-1+T-1, D12-1=D21-1=R-1-T-1, D22-1=R-1+T-1+Φ-1, d=12⊗R-1yq, 12´=11 and "⊗" denotes the Kronecker product between matrices (^{Langville y Stewart, 2004})..

Estimator of θ _q

From Equation 3, the empirical Bayesian estimator of θq´=aq´γq´ is:

(4)

The components of variance: σaq2, σγq2 and σeq2 can be estimated from the marginal distribution of y _q using restricted maximum likelihood (^{Lynch and Walsh, 1998}; ^{Vattikuti et al., 2012}).

The multivariate linear model

When two or more traits are used in the prediction of the breeding value, the combined linear model from Equation 1 can be written as:

(5)

where, now, y´=y1´ y2´… yt´ ~ NMV(0,V), a´=a1´ a2´…at´ ~ NMV(0,S), γ´=γ1´ γ2´…γt´ ~ NMV(0,Ω) y e´=e1´ e2´…et´ ~ NMV(0,Ψ) are vectors made up of t subvectors g×1 of observations (y), of additive genetic effects unexplained by the MM (a), of additive genomic effects (𝛾), and of errors (e), respectively; V = S + 3Ω + Ψ, where S = C ⊗ A, Ω = Γ ⊗ G and Ψ = E ⊗ I _g; C=σaqi (q,i = 1,2,…,t; t=number of traits) is the matrix of variances and covariances of the additive genetic effects unexplained by the MM (a), Γ=σ𝛾qi is the matrix of variances and covariances of additive genomic merits (𝛾), and E=σeqi is the matrix of variances and covariances of the residues; Z is an identity matrix (or of incidence) of order gt×gt; A, G and I _g are defined in Equation 1. The matrices C=σaqi and Γ=σγqi can be conformed from the estimations of the components of variances: σa2q, σγ2q and σe2q, and of the respective covariances (^{Vattikuti et al., 2012}).

Estimation of a and 𝛾

Let θ´=a´γ´, be a vector conformed by a´=a1´ a2´…at´ and γ´=γ1´ γ2´…γt´ (Equation 5); the posterior distribution of θ is similar to the distribution of θq´=aq´γq´ (Equation 3), thus, the empirical Bayesian estimator of θ is similar to Equation 4, i.e.,

(6)

where, now, the components that make up the D ^‒1 matrix are: D11-1=Ψ-1+(S-Ω)-1, D12-1=D21-1= Ψ-1-(S-Ω)-1 and D22-1=Ψ-1+(S-Ω)-1+Ω-1; d=12⊗Ψ-1y, Ψ-1=E-1⊗Ig and 12´=1 1.

Prediction of the breeding value in the first cycle of selection

In the first selection cycle the predictor of the breeding value of the candidates for selection (θ¯^) can be written as:

(7)

where â and ŷ are sub-vectors of θ^BE=Dd (Equation 6).

Prediction of the breeding value after the first selection cycle

In order to obtain the predicted values for the candidates for selection from the second selection cycle, it is necessary to estimate the values of the vector u´=u1´ u2´…ut´ in the training population from the equation γ = X _t u, where X _t = I _t ⊗ X, I _t is the identity matrix t×t and X is the matrix of coded MM values in the training population. An estimator of u in the training population is:

(8)

where ŷ is the sub-vector of Equation 6. From Equation 8, the empirical Bayesian predictor of the breeding value after the first selection cycle is:

(9)

where W _l = I _t ⊗ X _l (l=2,3,…,N; N=number of selection cycles), I _t was defined previously, and X _l is the matrix of coded MM values obtained in the selection cycle l. Thus, from the second selection cycle, the only thing that will change in Equation 9 will be the coded values in matrix X _l .

Criterion to compare the efficiency of the prediction models

Since accuracy is equal to the correlation between predicted and observed values, their maximum value is 1. Assume that ρ _c and ρ _g denote the accuracy in the combined and the genomic linear model, respectively, then:

(10)

is the efficiency (^{Bulmer, 1980}) of the combined linear model with regard to the genomic linear model. Thus, if p=0, the efficiency of both models is equal (ρ _c = ρ _g ); p>0 si ρ _c > ρ _g (the efficiency of the combined model is greater than that of the genomic model) and if ρ _c < ρ _g , p<0 (the efficiency of the combined model is lower than that of the genomic model). Thus, Equation 10 allows determining the most efficient linear model to predict the genetic merit.

The genomic model is nested in the combined model

One of the most important results in the theory of GS is that the expectation of the genomic relationships matrix G is equal to the numerical relationships matrix A, i.e. , E (G) = A (^{Habier et al., 2007}). This means that G is a particular realization of A and that when the number of MM and genotypes increases in the training population, the value of G will tend concentrated around A, so that it can be assumed that at the limit, G=A. The same is true with the additive genomic variances and covariance matrix Γ in relation to the additive genetic variances and covariances matrix C. That is, when the number of MM and genotypes increases, the matrix Γ approaches C, and at the limit, Γ=C. When G=A and Γ=C, S=Ω and the matrices that make up the matrix D-1:D11-1=Ψ-1+(S-Ω)-1, D12-1=D21-1=Ψ-1-(S-Ω)-1 and D22-1=Ψ-1+(S-Ω)-1+Ω-1, are reduced to D11-1=Ψ-1, D12-1=D21-1=Ψ-1 and D22-1=Ψ-1+Ω-1, and matrix D-1 will be equal to Ψ-1Ψ-1Ψ-1Ψ-1+Ω-1^‒1.

This indicates that all breeding value information is concentrated in the additive genomic effects 𝛾 and that the values of vector a are null. In such a case, the empirical Bayes estimator θ^BE=Dd (Equation 6) becomes the predictor of the additive genomic merit (ŷ) and can be denoted as:

(11)

This result indicates that the genomic linear model is a particular case of the combined linear model.

The model with only phenotypic information is nested in the combined model

When the information of the MM is not used, matrix Ω is null and, in such a case, θ^BE becomes the predictor of the additive genetic effects (â) and can be written as:

(12)

This shows that the linear model with only phenotypic information is a particular case of the combined linear model. The â will be called standard predictor.

Accuracy of the three prediction models

The predicted values of the breeding value of the candidates for selection associated to each one of the three traits of the two maize populations (Populations 1 and 2) and of the wheat population (Population 3) were denoted as θ¯^₁, θ¯^₂ and θ¯^₃, for the combined model (Equation 7); γ^1, γ^2 and γ^3 for the genomic model (Equation 11), and a^1, a^2 and a^3 for the standard model (Equation 12). With the predicted and the observed values the accuracy (correlation between the predicted and observed values) was calculated for each one of the three traits of the three models; these are shown in Table 1.

Table 1: Correlations obtained between the predicted values of the standard, genomic and combined models, and the values of the observations of three traits in two maize populations and one wheat population. 

♦ ^†Grain yield (Mg ha^‒1), ^❡Ear height (cm), ^§Plant height (cm).

Numerical evaluation of the three prediction models

The efficiency of the combined model with regard to the standard model and the genomic model; and the efficiency of the standard model with regard to the genomic model, was evaluated through Equation 10 with the correlation values presented in Table 1.

Maize population 1

Maize population 2

Population 3

Advantage of the genomic model with regard to the standard model

In the GS the usual way of predicting the plant and animal breeding values in breeding programs is to substitute the numerical relationships matrix (A) by the genomic relationships matrix (G) in the prediction equations. Therefore, the prediction equation of the genomic model (Equation 11) and the standard model (Equation 12), are formally equivalent. When the number of MM and genotypes is large, both models tend to provide predictions that are increasingly more similar (Table 1, Population 3). However, the advantage of the genomic model with regard to the standard model lies in the possibility of reducing the intervals between selection cycles in more than two thirds. Thus, the genomic model is more efficient than the standard model when the efficiency is measured per year and not per selection cycle. According to ^{Beyene et al. (2015)}, the genomic selection requires 1.5 years to complete a selection cycle, while the phenotypic selection requires 4 years for each selection cycle.

Importance of the combined model

There are several Bayesian (^{Gianola, 2013}) and non-Bayesian (^{VanRaden, 2008}) methods to predict the breeding value in the univariate context under the assumption that the number of genotypes and MM is sufficiently large in the base population. In practice, however, not all the candidates for selection (plants or animals) have molecular markers. Therefore, a model such as the one proposed could be easily adapted to this case, thus increasing the precision in the prediction.

Empirical Bayes compared with GBLUP

Because the MM effects have a multivariate normal distribution, empirical Bayes and GBLUP should give very similar results (^{Robinson, 1991}) when the same prediction model is used. This is because the assumptions of GBLUP and empirical Bayes are basically the same and because, when the variances of the parameters are known, GBLUP is considered a particular case of the Bayesian methods (^{Blasco, 2001}).

Finally, how could the breeding value be predicted? Through the empirical Bayes proposed, through GBLUP or with some of the existing Bayesian approximations? The standard Bayesian models provide a better control of the uncertainty associated to the prediction of the breeding value (^{de los Campos et al., 2013}, ^{Gianola, 2013}), which is attained with much computation work (^{Verbyla et al., 2009}). GBLUP, in turn, requires the knowledge of variances of the parameters; when these variances are unknown, the statistical properties of GBLUP are also unknown (Gianola, 2013). According to ^{Blasco (2001)}, the election of one prediction model over another one should be based on the fact that the chosen model offers a solution that the others do not, on how ease is to solve the problem, and on the trust in their results. This last point has the greatest importance, since if the researcher feels comfortable with a specific method, it means that he/she knows its limitations and advantages and knows what to expect from the model when using it in a specific statistical analysis.