Approximated prediction of genomic selection accuracy when reference and candidate populations are related

Genomic selection is still to be evaluated and optimized in many species. Mathematical modeling of selection schemes prior to their implementation is a classical and useful tool for that purpose. These models include formalization of a number of entities including the precision of the estimated breeding value. To model genomic selection schemes, equations that predict this reliability as a function of factors such as the size of the reference population, its diversity, its genetic distance from the group of selection candidates genotyped, number of markers and strength of linkage disequilibrium are needed. The present paper aims at exploring new approximations of this reliability. Two alternative approximations are proposed for the estimation of the reliability of genomic estimated breeding values (GEBV) in the case of non-independence between candidate and reference populations. Both were derived from the Taylor series heuristic approach suggested by Goddard in 2009. A numerical exploration of their properties showed that the series were not equivalent in terms of convergence to the exact reliability, that the approximations may overestimate the precision of GEBV and that they converged towards their theoretical expectations. Formulae derived for these approximations were simple to handle in the case of independent markers. A few parameters that describe the markers’ genotypic variability (allele frequencies, linkage disequilibrium) can be estimated from genomic data corresponding to the population of interest or after making assumptions about their distribution. When markers are not in linkage equilibrium, replacing the real number of markers and QTL by the “effective number of independent loci”, as proposed earlier is a practical solution. In this paper, we considered an alternative, i.e. an “equivalent number of independent loci” which would give a GEBV reliability for unrelated individuals by considering a sub-set of independent markers that is identical to the reliability obtained by considering the full set of markers. This paper is a further step towards the development of deterministic models that describe breeding plans based on the use of genomic information. Such deterministic models carry low computational burden, which allows design optimization through intensive numerical exploration.


Background
The effectiveness of genomic selection comes from the possibility of predicting breeding values on un-phenotyped and young animals [1]. Genomic selection promised and proved to be extremely efficient and beneficial for dairy cattle (e.g. [2][3][4][5][6][7]), but debate continues for other species and production sectors (e.g. [8][9][10][11][12]). A key criterion to decide whether or not selection schemes (also referred to here as breeding plans) should include genomic information is the reliability of the genomic predictor. It was clearly shown that this reliability depends on the structure of the reference population and on the characteristics of the marker set used. The size of this reference population, its diversity, the genetic distance between the reference and the group of selection candidates genotyped, the number of markers, and the degree or strength of the linkage disequilibrium are the main factors that influence this reliability [13][14][15][16][17][18][19][20][21][22][23].
An extensive literature exists on the mathematical modeling of selection schemes prior to their implementation, in order, for instance, to optimize their design, or to evaluate the usefulness of new technologies such as embryo transfer, sperm selection, DNA markers and others (e.g. [24] for a review). These models account for factors such as selection intensities and maintenance or loss of genetic variability. Among these parameters, the precision of breeding value estimates is central. To model genomic selection schemes, equations that predict this reliability as a function of the factors cited above are needed (e.g. [6,25,26]).
The quantitative influence of these factors (size of the reference population, its diversity, etc.) was assessed by simulation studies [18-21, 27, 28]. An equation that predicts the reliability of genomic evaluation in the very simple situation of independent quantitative trait loci (QTL), that are perfectly marked by single nucleotide polymorphisms (SNPs) and populations (reference and candidates) of unrelated individuals was derived [13]. This approach was extended to the case when only a part of the genetic variability is imperfectly marked by SNPs [16,19], and the situation of non-independence between reference and candidate populations was explored [17]. It was demonstrated that genomic information captures historical linkage disequilibrium, short-term linkage between QTL and markers and additive relationships between reference and candidate individuals, the equation of the reliability accounting for these three phenomena being derived in a very simple case of one QTL marked by a single SNP [22].
A Taylor expansion of a matrix inverse involved in the reliability formula was suggested [18], which led to the algebraic development of an approximation. This approximation seems to work well in the simple situation but lacks generality. In this paper, an alternative approximation is proposed, opening a way to include non-independence between reference and candidate populations, and between markers.
After a formalization of the genomic selection context, the principles that underlie these approximations are presented and their properties are compared by using a simple example. Then, the new approximation is derived when reference and candidate animals are related. This is illustrated by some numerical examples. Finally, the extension to the linkage disequilibrium situation is described.

General framework
Although the prediction equations derived below were based on a number of simplifying assumptions, it is important to first draw a complete description of the biological framework, as a basis to subsequently simplify the discussion.
The SNP effects are estimated in a reference population, P r , comprising n r individuals. The genomic estimated breeding values (GEBV) are calculated for a population of candidates for selection and used in breeding, P c , comprising n c individuals.
Let P = (P r , P c ) the population structure (including pedigree relationships between individuals and marker allele frequencies, but not including genotypes and phenotypes).
Individuals are characterized by their genotypes at n M markers (observed) and at n Q QTL (unknown). Alleles will be noted A m and B m for the marker m and A q and B q for the QTL q. Let a tim ∊ {0, 1, 2} and a tiq ∊ {0, 1, 2} be the numbers of B m (and respectively, B q ) alleles that an individual i from population P t (P r or P c ) carries at marker m (respectively, QTL q). Let p tm and p tq be the frequencies of alleles B m and B q in P t .
Genotypic values will be assigned to the different markers and QTL genotypes. Following [18], genotypes will be coded as x tim = a tim − 2p tm and w tiq = a tiq − 2p tq . Different codifications can be proposed [15]. In particular, as described for instance in [29], genotypic values may be standardized, i.e. x tim = (a tim − 2p tm )/σ tm and w tiq = (a tiq −2p tq )/σ tq , with variances σ tm 2 = 2p tm (1 − p tm ) and σ tq 2 = 2p tq (1 − p tq ). Most of the following developments are given with the first codification here, and the results with the second codification are described in a specific section.
These genotypic values are assembled in matrices X (dim (X) = (n r + n c ) × n M ) and W (dim (W) = (n r + n c ) × n Q ). Sub-matrices corresponding to sub-populations will be noted in the following way: . The genetic model assumes additivity of QTL effects. The additive genetic value of an individual is described as g ti = n Q q=1 w tiq α q and, in general, g = Wα. The phenotypic values when observed are y = g + ɛ.
A statistical model describes the performances in the reference population as random variables for which the expectations are linear combinations of SNP effects: y ri = n S m=1 X rim β m + e ri and, in general, y = Xβ + e. In these models, the SNP (or QTL) effects may be considered as fixed, or random. Since the number of SNPs is much bigger than the number of individuals (n M ≫ n r ) the second solution is generally chosen in the statistical model (but not always see [1,13]). In the random model, a distribution L θ β , V β (respectively L(θ α , V α )) of the SNP (respectively QTL) effects is assumed, with θ β (respectively, θ α ) being the vector of expectations and V β (respectively, V α ) being the matrix of variances. For a full description of the variability, the V β and V α matrices are each subdivided into four blocks corresponding to the reference and candidate populations and their covariances. Covariances between the α and β vectors have also to be considered. Most generally, the SNP (QTL) effects are supposed i.i.d. giving V β = Iσ β 2 (V α = Iσ α 2 ). The interpretation of these QTL effects is nicely debated in Gianola et al. [30]. In the frequentist view, we simply have to imagine that QTL effects are randomly sampled from a distribution with a σ α 2 variance. In the Bayesian context, the prior variability of the SNP effects was most generally described as heteroskedastic or even coming from mixtures of SNPs with or without an effect on the trait.
The expectations θ β (θ α ) are generally assumed equal to zero, but when information about population history is available (in particular, when we know it is a mixed population), non-zero values should be considered.
The vector q = Xβ is a quantity similar but not equal to the genetic value g. Its element q ti is the molecular score of individual i in population t. This vector may be segmented in two parts: q ′ = q ′ r , q ′ c . Since the variances may be defined within a population, Assuming that the distribution of marker effects is centered (θ β = 0) and i.i.d. (V βr = Iσ 2 βr and V βc = Iσ 2 βc ), and extending Gianola et al. [30] βr m σ 2 mr = σ 2 βr τ r in the reference population, and v(q ci ) = σ 2 βc τ c in the candidate population. Assuming that the distribution of the marker effects and genotypes are the same in P r and P c , i.e. p rm = p cm = p m , p rq = p cq = p q , thus τ r = τ C = τ and σ βr 2 = σ βc 2 = σ β 2 , we define σ q 2 = τσ β 2 . Thus, v(q|X) = 1 τ XX ′ σ 2 q . These equations hold even if the markers are in linkage disequilibrium (LD) as shown in Eq. A2 from Gianola et al. [30].
We note σ 2 as the total phenotypic variance, i.e. σ 2 = σ q 2 + σ e 2 , and ν 2 as the proportion of this variance explained by the molecular score will be noted γ.
The SNP effects β may be estimated in different ways. The genomic best linear unbiased prediction (BLUP) will only be considered here, with β = cov β, y var y −1 y .
Classically, this equation becomes The linear combination q c = X cβ is the GBLUP vector for candidates in P c . It must be emphasized that these estimations and predictions are conditional on the genotypic structures defined by X (X r and X c ).
Given X, the reliability of the GBLUP is In [16], the reliability is described (Eq. 6 in [16]) as r g ci ,q ci = r g ci , q ci × r q ci ,q ci , by ignoring the conditioning on X. In Goddard et al. [18], the reliability is described as is the proportion of the genetic variance explained by the markers and v(q ci) v(q ci ) is the accuracy of estimated marker effects. This is similar to the qrQ reported by Dekkers et al. [25]. All these reliability formulae are approximations since cov 2 g ci ,q ci = cov 2 w ciq α q , x cisβs � = v q ci = v x cisβs , in general.

Situation analyzed in this paper
In the following, ignoring the difficulty that was mentioned above, we will assume r 2 q ci , v(q ci |X) . We are interested in a single candidate in P c with a x c vector of marker genotypes.
Formulae were simplified in two ways. (1) the i index of the candidate was omitted in the following developments: the genetic value of the candidate is noted q c , estimated by q c = cov q c , y v y −1 y, and its precision where x c is a row vector); (2) the r indices of reference individuals were most often omitted, which resulted in y i for their phenotypes and q i for their molecular scores.
In fact, our objective was to estimate the expectation of this precision across the variation domain of X r and x c given the pedigree structure (P r , P c ) : E X r 2 q c ,q c |X |P . It will be noted E r 2 q c ,q c . The following approximation was made: . Let A be the pedigree relationship matrix between individuals in P. Its blocks are A = a cc A cr A rc A rr .
It must be noted that the σ e 2 term in the diagonal of the V submatrix corresponding to the candidate population is artificial since candidates are not phenotyped.
We have E[G*] = Aτ. The limits of this equality will be discussed below. As indicated above, the denominator of the expected reliability E X [v(q c |X)], is τσ β 2 = σ q

2
. Approximating E v q c by E[cov(q c , y)]E[v(y)] −1 E[cov(y, q c )] is useless because it makes an oversimplification of the relationships between the reference and the candidate population: it considers separately the marginal distributions of x c X r ′ and (X r X r ′ + Iλ β ) −1 , while these random matrices are correlated. Estimating directly E[cov(q c , y)v(y) −1 cov(y, q c )] seems impossible in the general case. The approach of Goddard et al. [18] avoids this difficulty, i.e.
is approximated by a second degree Taylor expansion

Alternative approximations of the reliability Extension of Goddard's formula
In their "heuristic approximation for V *−1 ", Goddard et al. [18] considered the situation where unrelated individuals are included in the reference and candidate populations, that is E[G * ] = Iτ and G * = Iτ + E, with E, a "noise" matrix centered on the null matrix 0. A direct extension of their development would be the following. The matrix V = 1 τ G * σ 2 q + Iσ 2 e can be written as: The inverse matrix [I + Dγ] −1 will be approximated using a Taylor series. It must be emphasized that the Taylor series I − Dγ + (Dγ) 2 − (Dγ) 3 Finally, the reliability of the candidate GBLUP is approximated by: A difficulty with this approximation comes from the T term. As an example, consider a reference population composed of n r half-sibs of the candidate, T = ξI + ψJ with ξ = 4 4+3γ . As T t = ξ t I + n t r ξ t + · · · J, the J coefficient will tend to ∞ as soon as n r ξ = 4n r 4+3γ > 1, a very realistic situation. Thus, the convergence of the Taylor series will be a balance between the increase of T t and decrease of [Dγ ] t . (1)

Principle
Using the classical matrix inversion lemma, the variance v q c |x c , X r = σ 2 β x c X ′ r X r X ′ r + I β −1 X r x ′ c may also be defined as v q c |x c , r X r is a very large matrix (n M × n M ) that describes the LD between markers: its elements tend to be smaller when they are more distant from the diagonal.
Elements of E[X r ′ X r ] are the following: E X ′ r X r ml = E i (a im − 2p m )(a il − 2p l ) = 2n r � ml , with ml the LD between loci m and l.
. . , σ 2 n M , the X r ′ X r + Iλ β matrix may be written as: which results in: Thus, we expect the Taylor series to converge to

First order approximation
c and the expectation of the reliability of the candidate GBLUP is approximated by Using Finally, the expectation is: corresponding exponents are summed (e.g. α iiil 1111 = α il 31 ). The resulting X im moments are in Table 2.

Second approximation
The expectations E[x 2 cm x 2 rim ] and E x 2 cm x rim x cl x ril are also obtained from the coefficients in Table 1, i.e.: ci and, when markers . After some algebra, it appears that: where ā 2 ci , ᾱ 22 ci and γ 22 ci are the means of the corresponding coefficients, considering all possible i reference individuals.

Parameter estimation
The from Additional file 1 δ s are the 15 classical identity states probabilities between two individuals [33][34][35] Coefficients of expectations E X i X j X 2 k and E[X i X j X k X l ] involve IBD status between three or four different individuals and are explained in Additional file 1

Application in the case of independent markers
This situation either assumes low density marker information, or corresponds to the idea of an effective number of loci that was developed by Goddard [16,31]. In the first case, the proportion of the genetic variance explained by the markers v(q ci )

v(g ci)
is small and this quantity should be considered when estimating the genomic precision.

First approximation
Using the notation a ij the coancestry coefficient between individuals i and j, are functions of the probabilities of the identity states between gametes of ij · · · K individuals at marker m (Table 1). In the summations above, when individuals are repeated (e.g. i = j), the unknown. Their expectations can be derived by making assumptions about the distribution of the marker allele frequencies. They were derived assuming either a uniform distribution of allele frequencies or the U-shaped distribution of allelic frequencies proposed by Goddard [16]: with the constant k estimated as 1/ log2N e , N e being the effective size of the reference population. The expectations of the parameters are in Table 3. The corresponding algebra is detailed in Additional file 2.
The parameters τ and τ 2 are linked to the number M e of independent segments. This quantity M e was defined by Goddard [16] as the number of independent chromosomal segments which would give the same variance of genomic covariances c ij between individuals i and j as that observed, i.e. when LD exists. Conditional on the genotypic observation, the genomic covariance between two individuals is cov( . When the markers are in linkage equilibrium, the covariance term is null, From the appendix in the paper of Goddard [16], this variance is v X (c ij ) = σ q 4 /M e . Thus: It must be emphasized that M e , which depends on the variability of allele frequencies, is not the number of markers n M .

The case of unrelated individuals
The first approximation gives results similar to Goddard et al. [18] when individuals are not related. In this case, is the proportion of the phenotypic variance explained by the molecular score.

Non-independence between reference and candidate population, a simple example
We consider the situation of a candidate that is the son of one of the n r individuals in P r (say the first in the list) while still assuming that reference individuals are unrelated. In this situation, the pedigree relationship matrix is . Applications of formulae (2) and (3) are described in Additional file 3. The expected approximate precision with the first approach is: where c1 = (a + b) 3 + a 2 + b 2 1 2 a, c2 = 1 4 a b 2 − a 2 and c3 = a 2 + b 2 + 1 2 ab. And with the second approach:

Alternative genotypes codification
In all the previous developments, genotypes were coded x tim = a tim − 2p tm and w tiq = a tiq − 2p tq . Alternatively, we could define x tim = (a tim − 2p tm )/σ tm and w tiq = (a tiq − 2p tq )/σ tq . The relation between genetic and marker variances becomes σ q 2 = n M σ β 2 and the relation between pedigree and genomic matrices becomes E[G*] = An M . Thus, formulae (1) and (2)

The case of markers in linkage disequilibrium
So far, following Goddard [16], we considered the situation of n M independent segments that each carries a single QTL in LD with a single marker. More typically, the genomic information consists of a large number of nonindependent markers. This non-independence comes from long-term effects due to bottlenecks, mutations, migrations, etc. and short-term effects due to family structure.

Effective and equivalent numbers of independent loci
We based our developments on the very fruitful concept of the effective number of loci that Goddard defined as "the number of independent loci that gives the same variance of realized relationships as that obtained in the more realistic situation" (Goddard [16] appendix). Since our objective was to predict the reliability of GEBV, we now suggest the alternative definition of an "equivalent number of independent loci" which would give the reliability of GEBV for unrelated individuals when considering a sub-set of independent markers that would be identical to the reliability obtained when considering the full set of markers. From the derivation of the reliability given previously, defining x c i and X i r as the genotype matrices of the independent loci, we need With a few simplifying assumptions (identical distribution of genotypes in the reference and candidate populations and equal genotypic variance at all loci) a simple formula can be derived [see Additional file 4]: where tr[M] is the trace of matrix M.
Once marker allele frequencies and between-marker LD are estimated in a population of interest, the equivalent number of independent loci which can be estimated from formula (9) and this parameter can be used in models that predict the genetic gain expected from a genomic selection scheme applied to this population.
In the more general situation, prior to the observation of the X r matrix, a simple approximation for n M i is obtained assuming equal variances σ m 2 = s 2 , and using the relation between expected LD and effective population size N e as derived by Sved [32]: .

Towards an exact treatment of linkage disequilibrium
For a complete treatment of the LD situation, it is necessary to estimate the expectations of the product of four genetic values. For instance, with the second approximation [formula (2)], we need to compute E[x cm X rim x cl X ril ]. Let X im = g imf + g imd , where g imf and g imd are the "values" of the alleles transmitted to individual i by its father and its dam, with g imf and g imd = (0 or 1) − p m . They will be called allelic values in the following. Equivalent terms are defined for x cl ,x cm and X il . The random variable M cls is the allele of individual c received from s at locus l f or d . M cmt , M ilu and M imv are defined similarly.
For the candidate c as for the reference individual i, the pair of genetic values may originate from the same parent (and coded on the same chromosome) or not, giving four types of g cls , g cmt , g ilu, g imv vectors. In type 1 (s = t and u = v), both alleles (belonging to loci m and l ) of each pair of loci (one for c and one for i) are on the same chromosome (may be from the two fathers, the two dams, c's father and i's dam or i's father and c's dam). In type 2 (s = t and u � = v), both alleles (belonging to loci m and l) of the candidate are on the same chromosome, while alleles of the reference individual i are not on the same chromosome.Type 3 (s � = t and u = v) is the reverse from type 2. In type 4 (s � = t and u � = v), alleles of loci m and l of both individuals and i are on different chromosomes.
For each of these situations, the identity by descent (IBD) status between alleles at locus m on chromosomes ct and iv, and at locus l on chromosomes cs and iu are considered. There are four, as follows: Thus, 16 terms involved in E[x cl x cm X il X im ] are given by: As described in Additional file 5, only seven E[g cls g cmt g ilu g imv |S k ] are non-null (Table 4). Principles on which the probabilities φ k stuv are estimated and basic examples are described in Additional file 5.
As an illustration, we consider again the situation of a candidate (c), that is the son of one of the n r individuals in P r and assume that c 's dam is unrelated to the sire. In formula (2), the summation over the reference individuals i comprises a single term for the sire of the candidate and n r − 1 terms for the individual that are unrelated to the c members of this reference population.
Based on Additional files 1 and 5, expectations involved in the precision formulae (2) are:

Basic situation: no LD and unrelated individuals
Convergence of Taylor series and quality of expectation of the reliability approximations were tested for different population sizes (n r = 500, 1000, 1500 and 2500), numbers of markers (n M = 50, 100, 250, 1000, 1500, 2000 and 2500) and proportions of the phenotypic variance explained by the molecular score (ν 2 = 0.1, 0.4 and 0.7). Given the set of allele frequencies p m (m = 1 . . . n M ), genotypes X of n r + 1 individuals were generated and the G matrix was built. The reliability of the candidate individual GEBV, r 2 = v(q c |X) v(q c |X) was computed as described in the section «Situation analyzed» , as well as approximations considering 1-10 elements in the Taylor series I − Dγ + D 2 γ 2 − D 3 γ 3 ··· The convergence of the series as predicted by the value (lower or higher than 1) of the matrix's largest eigenvalue was checked numerically, by estimating the mean of this largest eigenvalue from five simulations in each case studied (n r = 200 to 1000; n M = 100 to 2000 and ν 2 = 0.1, 0.4, 0.7).
This limited number of replications was chosen after observation of a very limited variance of this eigenvalue. Finally, the asymptotic values of the suggested approximations [formulae (5) and (6)] were computed using the number of independent segments as described by [4]. The process was repeated 50 times and the means of those exact or approximated reliabilities computed. Figure 1a and b illustrates the convergence of the Taylor series when 2000 markers are used, and Tables 5 and  6 give the results for both approximations when ν 2 = 0.4. The Taylor series converged when the proportion ν 2 of the phenotypic variance explained by the molecular score was low, with oscillations and divergence observed when ν 2 = 0.4 or 0.7 with the first approximation and ν 2 = 0.7 with the second approximation. These observations were in accordance with the deviation to one of the largest eigenvalue of the matrix involved in the series (Fig. 2a, b). When the series converged, the approximations rapidly reached a plateau, at the 3rd (respectively, 2nd) order for the first (respectively, second) approximation. Table 6 shows that the second Taylor series converges always when ν 2 = 0.4. The proposed approximation was generally biased upwards. This over-estimation of the precision was generally limited but increased as the number of markers and the reference population size decreased. The maximum over-estimation observed was

Table 4 Expectations of products of four allelic values received by two individuals at two loci depending on the IBD status and parental origins of the alleles
Only non-null terms are given p m and p l are the frequencies of the most frequent alleles at loci m and l. lm is the linkage disequilibrium measure between m and l g cls = (0 or 1) − p l is the allelic value the candidate received from its parent s at locus l etc S ml means c and i genes are IBD at m and l, only at m etc 37.5 % (0.22 instead of 0.16 with a standard error less than 0.02). Based on the results in Table 5, it appears that when the first Taylor series converges, the proposed approximation is also slightly over-estimated. The expectation of the approximations, as given in formulae (5) and (6) are very close to the observations.

No LD and non-independence between reference and candidate population
The quality of the approximation was tested as above, by considering the case of a candidate having one of its parents in the reference population and all other individuals being unrelated. Tables 7 and 8, which summarize the results of the simulation, show that the second approximation is still the most efficient (systematic convergence of the Taylor series and consistency between first order approximation and its expectation). Again, an overestimation of about 20 % is observed with this approximation.

Example of the use of the second approach
As an illustration of formula (3) different situations that differ in the relationships between the candidate and reference populations were compared. Coefficients of formula (3) were estimated using the elements in Table 3. An effective reference population size of 200, the genotyping of 10,000 markers and a heritability of 0.4 were assumed. Scenarios included no individuals related to the candidate in the reference population, its sire, both parents, 1-10 half-sibs (or uncles), and a combination of parental and half-sib information.
The results are in Fig. 3. The linearity of the precision increases with the number of half-sibs, which is consistent with the approximation, but unsatisfactory, as discussed below.

Equivalent number of independent loci
This number was computed using formula (8), for various effective population sizes (N e = 100 to 1000 ), heritabilities (h 2 = 0.1 to 0.5), total numbers of loci Table 5 Performances of the first approximation r 2 q c ,q c for an unrelated reference population as a function of the number of markers (n M ) and reference population size (n R ), assuming ν 2 = 0.4 The convergence criterion is the value of the Taylor series at order 10 E r 2 qc,qc is the expectation of the first approximation across the distribution of allele frequencies as given in Goddard [16]  (n M = 1000 to 10, 000) and reference population sizes (n R = 1000 to 2500). Figure 4 shows how equivalent numbers of independent loci (n M i ) vary with the total number of markers (n M ) and reference population size (n R ). As n M increases, the number n M i rapidly converges to a value which strongly depends on the size of the reference population size. This dependence on n R of the equivalent number of independent loci does not exist in the Goddard's effective number of loci and clearly shows the difference in nature between these concepts. Three phenomena, observed when considering the extreme case of two markers (see Additional file 5), explain this behavior: (1) The trace T of (E[X r ′ X r ] + λ β I) −1 is a decreasing function of n r : as a consequence, the larger is the population size, the smaller is T, which is proportional to the marker effects conditional variances v(β) − cov β, y v(y) −1 cov(y, β)) and the higher is the variance of the estimated molecular score (v q c |y = x c cov β, y v(y) −1 cov y, β x ′ c ).
(2) The trace T is always higher in the situation of LD than for independent markers (T LD > T LE ). (3) The rate of decrease is higher for T LD than for T LE .
On the whole, the reliability for a given number of observed markers corresponds to the reliability that is reached with a larger number of independent loci when the size of the reference population is larger. Figure 5 indicates that the equivalent number of independent loci increases with heritability and effective population size. This last observation was expected since with larger effective population sizes, the LD between two loci decreases and this increases the effective number of loci. The effect of heritability is less direct.

Discussion
The objective of this paper was to explore approximations of the precision of genomic selection when the selection candidate has relatives in the reference population. Two approximations were developed and numerically compared.
These approximations were based on Taylor expansions of a matrix inverse M −1 . In both cases, the initial matrix is the sum of the identity matrix and a perturbation (M = I + E). Convergence of these series is not guaranteed and depends on the behavior of the perturbation With the first approximation, derived from the appendix in [18], this convergence failed when the number of makers was too small (less than 1500 in our example) or the heritability was greater than 0.1. This was only observed when ν 2 = 0.7 with the second approximation. This is fully consistent with the deviation to one of the largest eigenvalue of the E matrix.
The expectation of the proposed approximation, when data were simulated with the model corresponding to the hypotheses underlying their algebraic development, was very close to the mean value after 50 simulations. Thus, extremely fast estimation of the precision is possible, which allows intensive optimization and comparison of selection schemes.
When individuals are unrelated and markers are in linkage equilibrium, we obtain an estimation of the GEBV accuracy which differs from that of Goddard et al. [18]. Table 6 Performances of the second approximation r 2 q c ,q c for an unrelated reference population as a function of the number of markers (n M ) and reference population size (n R ), assuming ν 2 = 0.4 The convergence criterion is the value of the Taylor series at order 10 E r 2 qc,qc is the expectation of the second approximation across the distribution of allele frequencies as given in Goddard [16]  This is surprising since that approach was said to be based on the Taylor approximation used here. Their formula may be obtained in a simpler way [see Additional file 6]. However, relaxing the assumption of "absence of between-individual relationship" is not straightforward using this approach. A strong limit of our new approximation comes from the limitation to the first order term of the Taylor series. Deriving algebra was only possible at this stage. The side effect is that no genotypic covariance terms between reference individuals appear in this approximation. As a consequence, only the direct relationships between candidate and reference individuals play a role in the estimation, but not the structure within the reference population. This is unfortunate, because accuracies of genomic prediction are obviously affected by the construction of the reference population. Our last numerical example, in which there is a linear trend with the number of half-sibs, reveals this drawback: two half-sibs of the candidates are treated as unrelated and the information that they carry is just the double of that of a single half-sib. Future developments should focus on this limitation, for instance to derive the expectation of the x c C −1 B 2 x c ′ term. The U-shaped density function f (p) of allele frequencies was defined as in [16] . A Beta distribution B(φ a , φ b ) for the allele frequencies was assumed by Gianola et al. [30], following Wright [34]. Assuming that the frequency distribution is centered on 0.5, i.e. Φ a = Φ b = Φ, this quantity Φ can be adjusted to fit the distribution of Goddard. Using the Chi 2 test as a fitting option, we observed that the adjusted φ rapidly decreased as the population size increased (Fig. 6), with a slower and slower evolution as the population size grew larger (with n r = 200, 000 the adjusted φ is 0.9750000). Using a Beta distribution could give more generality to the results. If the expectation of τ and τ 2 are easily derived from the moments generating function of Beta distribution (E[τ ] = n r a 2a+1 Table 7 Performances of the first approximation r 2 q c ,q c when the parents of candidate belong to the reference population as a function of the number of markers (n M ) and reference population size (n R ), assuming ν 2 = 0.4 The convergence criterion is the value of the Taylor series at order 10 E r 2 qc,qc is the expectation of the first approximation across the distribution of allele frequencies as given in Goddard [16] Table 8 Performances of the second approximation r 2 q c ,q c when the parents of the candidates belong to the reference population as a function of the number of markers (n M ) and reference population size (n R ), assuming ν 2 = 0.4 The convergence criterion is the value of the Taylor series at order 10 E r 2 qc,qc is the expectation of the second approximation across the distribution of allele frequencies as given in Goddard [16]   is not simple. However, these quantities are quite easily obtained by numerical integration. Thus, adjusting a Beta distribution to observed allele frequencies and numerically computing formula (3) parameters would be a feasible and more versatile implementation of our second genomic precision approximation.
Our work focused on the BLUP precision of the molecular score r 2 q ci ,q ci |X = v(q ci) v(q ci ) but left aside the proportion of the genetic variance that is captured by the markers v(q ci ) v(g ci ) . This last term could be treated as in Goddard et al. [18]: v(q ci ) v(g ci ) = b = n M n M +M e with M e the number of independent segments. As noted in the section on the general framework, the quantity v(q ci ) v(g ci ) = b × r(q ci ,q ci |X) is only an approximation of these GEBV reliabilities i.e. r 2 g ci ,q ci |X = cov 2 (gci,qci|X) v(g ci |X)v(q ci |X) . Equality between those quantities is obtained when X = W(identity between statistical and genetical models), a condition assumed in Goddard [16] where markers and QTL are modeled as a series of uncorrelated pairs.
All the developments shown in this paper are based on the hypothesis that the reliability of GEBV based on non-independent markers for a trait controlled by n Q QTL that are in incomplete LD with the markers can be approached by the reliability of GEBV when there are n M independent segments that carry a single QTL in LD with a single marker. A few difficulties arose when applying this approach proposed by Goddard [16]. How many independent markers should be considered? The reasoning in Goddard [16] was based on the idea of an effective number of loci (M e ) corresponding to a given variance of realized relationships. Here, we proposed the alternative equivalent number of independent loci (M i ) which corresponds to a given reliability. We showed that this M i number depends on the size of the reference population and on heritability, a dependence that does not occur with M e . If we invert the argument, controlling the level of realized relationships variance with the effective number of loci (M e ) does not seem to be a good approach to control the estimated GEBV reliability.
As detailed by Hayes et al. [17], the effective number of independent chromosome segments depends on the population structure. The higher is the mean relationship level, the smaller is this effective number. However, we suggest the use of this number as estimated from a set of unrelated individuals, or of its expectation prior to any observation, assuming independence between individuals. Without formal proof, the idea was that long-term LD was considered by using an effective (or equivalent) number of independent loci while short-term non-independence was taken into account with our formalization of the matrix's expectations that is developed in Additional file 1. A complete proof of the procedure is still needed.
Regardless of the definition of M e or M i , there is no reason that the number of independent loci must equal the number of QTL, which is unknown, contrary to the hypothesis about pairs of marker-QTL (in practice, since the QTL effects are random variables, many segments will only have very small effects on the trait, thus simulating the more likely situation of a limited number of "real" QTL). Equating X and W as well as σ β 2 and σ α 2 has no clear justification. The variance v q c |X of the molecular score should not be σ β 2 x c X r ′ (X r X r ′ + Iλ β ) −1 X r x c ′ but This other formula assembles two sets of unknown parameters: the variances σ α 2 and σ β 2 , and the genotypes X and W. It is often assumed that σ 2 β = σ 2 g /(n Mτ ) (e.g. [1]), which results in an overestimation of the β parameter since LD is not considered. Working on the number of independent loci (M e or M i ) apparently solves this difficulty. The QTL variance σ 2 α = σ 2 g / n Qτ could be derived based on a hypothesis about the number of QTL. The situation is more difficult for the genotype matrices since the W r matrix is not observed.
If the framework considered so far (n M markers-QTL pairs with strong LD within pairs and no LD between pairs) is partly retained, a slight improvement is possible considering the element b of the genetic variability  Fig. 6 Parameter of the beta distribution B(φ, φ) that best fits Godard's distribution of allele frequencies explained by SNPs. The idea would be to replace, in the formulae used in this paper, σ q 2 by b × σ g 2 . Element b can be derived by considering that the markers' (β) and QTL' (α) effects are fixed in the genetic and statistical models. Leaving aside the singularity of X r ′ X r when the number of SNPs is large, the marker effects are now estimated by β = X ′ r X r −1 X ′ r y and the molecular score defined as q = X rβ , while the genetic value was g = W r α. Given the genotype matrices, the sample genetic variability is v g = α ′ W r ′ W r α and the sample molecular score variability y ′ X r (X r ′ X r ) −1 X r ′ y with an expectation v q = α ′ W r ′ X r (X r ′ X r ) −1 X r ′ W r α. The part of the genetic variability explained by the SNPs is the ratio b = v q /v g .
Expectations of the matrices' product elements X ′ r X r ml are 2n r ml off diagonal and 2n r p m (1 − p m ) = n r σ m 2 in the diagonal, with similar expressions for W ′ r X r and W r ′ W r elements. Following Goddard [16], approximating expectations of the matrices' functions by the function of their expectation, and assuming that (1) markers are independent, (2) each QTL q is in LD with only one marker m(q), with a LD value � qm(q) , and (3) individuals are unrelated: v g = n r ∑ q α q 2 σ q 2 , we get v q ∼ 4n r q � 2 qm(q) σ 2 m(q) α 2 q = n r q r 2 qm(q) α 2 q σ 2 q , and b = q r 2 qm(q) α 2 q σ 2 q q α 2 q σ 2 q , corresponding to Eq. (4) in [16].
The ratio b is the weighted mean of LD r 2 . Unfortunately, neither α q 2 nor σ q 2 are known. The unweighted mean q r 2 qm(q) n q =r 2 may be a fruitful approximation. Following Sved [33], the expectation of r 2 qm(q) is 1 1+4N e c with c being the distance, in Morgan, between the QTL and its marker. Let L be the total length of the genome, and assume an equal distance L/n M between each successive marker b ∼ ) could not be used here due to algebra complexity. However, in the case of unrelated individuals and independent markers, numerical evaluation of the difference between exact and approximated results for various reference population sizes and numbers of markers shows a very small underestimation of the reliability ( Table 9).
The theory presented here was developed by considering a single selection candidate. When candidates are diversely related to the reference population, as suggested in Goddard et al. [18], the candidates should be examined one by one. Moreover, non-independence between candidates should be considered. A further step towards the modeling of genomic selection could be an approximation of the mean genetic values of selected individuals when GEBV reliabilities are heterogeneous.
A few other hypotheses were made in this paper, including additivity and i.i.d. of QTL effects, and the use of GBLUP. As long as the objective is to model and optimize breeding plans, only relative values are interesting and we assumed that these hypotheses were not critical.

Conclusions
The objective of this paper was to provide a further step towards the development of deterministic models that describe genomic breeding plans. Such deterministic models carry low computational burden and thus allow design optimization through intensive numerical exploration.
We proposed two alternative approximations of the estimation of GEBV reliability in the case of nonindependence between candidate and reference populations. Both were derived from the Taylor series heuristic approach suggested by Goddard [16]. A numerical exploration of their properties showed that the series were not equivalent in terms of convergence to the exact reliability, that the approximations may overestimate GEBV precision and that they perfectly converged toward their theoretical expectations. Formulae derived for these approximations were simple to handle in the case of independent markers. A few parameters that describe the markers' genotypic variability (allele frequencies, linkage disequilibrium) can be estimated from genomic data corresponding to the population of interest or estimated after assumption about their distribution.
When markers are not in linkage equilibrium (i.e. there is LD), replacing the real number of markers and QTL by an effective or equivalent number of independent loci, as proposed by Goddard [16] and Hayes et al. [17], is a practical solution. Research efforts are still needed to overcome some strong limits of this approach.
(métaprogramme Selgen). Andrew Swan, Julius van der Werf, Mike Goddard, Anne Ricard and Bruno Goffinet are thanked for their many useful comments. Rob Banks is particularly thanked for his help at many levels.