Metafounders are related to F st fixation indices and reduce bias in single-step genomic evaluations.

BACKGROUND: Metafounders are pseudo-individuals that encapsulate genetic heterozygosity and relationships within and across base pedigree populations, i.e. ancestral populations. This work addresses the estimation and usefulness of metafounder relationships in single-step genomic best linear unbiased prediction (ssGBLUP).
RESULTS: We show that ancestral relationship parameters are proportional to standardized covariances of base allelic frequencies across populations, such as [Formula: see text] fixation indexes. These covariances of base allelic frequencies can be estimated from marker genotypes of related recent individuals and pedigree. Simple methods for their estimation include naïve computation of allele frequencies from marker genotypes or a method of moments that equates average pedigree-based and marker-based relationships. Complex methods include generalized least squares (best linear unbiased estimator (BLUE)) or maximum likelihood based on pedigree relationships. To our knowledge, methods to infer [Formula: see text] coefficients from marker data have not been developed for related individuals. We derived a genomic relationship matrix, compatible with pedigree relationships, that is constructed as a cross-product of {-1,0,1} codes and that is equivalent (apart from scale factors) to an identity-by-state relationship matrix at genome-wide markers. Using a simulation with a single population under selection in which only males and youngest animals are genotyped, we observed that generalized least squares or maximum likelihood gave accurate and unbiased estimates of the ancestral relationship parameter (true value: 0.40) whereas the naïve method and the method of moments were biased (average estimates of 0.43 and 0.35). We also observed that genomic evaluation by ssGBLUP using metafounders was less biased in terms of estimates of genetic trend (bias of 0.01 instead of 0.12), resulted in less overdispersed (0.94 instead of 0.99) and as accurate (0.74) estimates of breeding values than ssGBLUP without metafounders and provided consistent estimates of heritability.
CONCLUSIONS: Estimation of metafounder relationships can be achieved using BLUP-like methods with pedigree and markers. Inclusion of metafounder relationships reduces bias of genomic predictions with no loss in accuracy.

Background

Metafounders are pseudo-individuals that describe relationships within and across pedigree base populations. The concept of metafounders provides a coherent framework for the theory of genomic evaluation [1]. Genomic evaluation in agricultural species often implies partially genotyped populations, i.e. some individuals are genotyped, using high-density genetic markers across the genome, others are not, and phenotypes may be recorded in either of the two subsets. An integrated solution called single-step genomic best linear unbiased prediction (ssGBLUP) has been proposed [2-4]. This solution uses the following integrated relationship matrix:with inverse:where is the genomic relationship matrix, is the pedigree-based relationship matrix, and matrices are submatrices of with labels 1 and 2 denoting non-genotyped and genotyped individuals, respectively.

Because genotyped animals are often not a random sample from the analyzed populations (they tend to be younger or selected), it was quickly acknowledged that a proper analysis requires specifying different means for genotyped and non-genotyped individuals for the trait under consideration. These different means can be considered as parameters of the model, which are either fixed [4] or random [5, 6] effects. In the latter case, the random variables induce covariances between individuals, a situation that is informally referred to as “compatibility” of genomic and pedigree relationships. In fact, compatibility implies equality of the average breeding value of the base population and of the genetic variance [7] across the different measures of relationships. Numerically, the problem appears as follows. The formulae for matrix and its inverse contain () and () (assuming is full rank), respectively. This suggests that if and differ too much, biases may appear.

Genomic relationships are usually computed in one of two manners: using “cross-products” [8] or “corrected identity-by-state (IBS)” [9]. Both depend critically on assumed allele frequencies at markers in the pedigree base population [10]. Base allele frequencies are often unavailable. However, for most purposes, allele frequencies are not of interest per se and can be treated as nuisance parameters that can be marginalized. Christensen [11] achieved an algebraic integration of allele frequencies, leading to a very simple covariance structure with allele frequencies in genomic relationships fixed at 0.5 (e.g., using genotypes coded as {−1,0,1} in the cross-product method) and a parameter called that describes relationships between pedigree founders i.e. in the base population. A second parameter in Christensen’s marginalisation is , which is a measure of marker heterozygosity in the base population. Therefore, instead of inferring (thousands of) base allele frequencies, inference can be based on two simple parameters and . Both can be estimated by maximizing the likelihood of observed genotypes. In addition, this approach considers the fact that pedigree depth is arbitrary and mostly based on historical availability of records.

Legarra et al. [1] showed the equivalence of the Christensen approach to the metafounder concept: pseudo-individuals that encapsulate three ideas: (a) separate means for each base population [4, 12, 13], (b) randomness of these means [5] and (c) propagation of the randomness of these means to the progeny [11], while accommodating several populations with complex crosses, e.g. [14]. Legarra et al. [1] also generalized one relationship between founders (scalar ) to several relationships between founders in the pedigree, i.e. ancestral relationships (matrix ), and suggested simple methods to estimate them. Legarra et al. [1] showed that construction of from and a pedigree reduces to the use of the tabular rules [15] for construction of relationships, and its inversion is achieved by inversion of followed by Henderson’s rules [16]. We provide an example of matrices and in “Appendix”. However, the performance of their model has not been tested so far, either for estimation of ancestral relationships or for genomic evaluation.

This work has two objectives. The first is to show that the structure of the metafounder approach yields an alternative parameterization and method for estimation of ancestral relationships. By doing so, we found that ancestral relationships are generalizations of Wright’s fixation index [17]. The second goal is to test, by simulation, (1) methods to estimate ancestral relationship parameters, (2) the quality of genomic predictions using metafounders, and (3) the quality of variance component estimation. For the second goal, the simulated population is undergoing selection and with a complete partially genotyped pedigree.

Methods

Relationship between metafounders and allele frequencies in the pedigree base population

Single population

Let be a matrix of genotypes coded as gene content, i.e. {0,1,2} and the genomic relationship matrix with a matrix of 1s, with reference alleles taken at random, so that the expected allele frequency is 0.5 for a random locus [11]. In other words, the matrix contains values of {−1,0,1} for each genotype. In a single population, let be the relationship coefficient between pedigree founders or, equivalently, the self-relationship of the metafounder [1, 11]. Parameter (defined above) is a measure of marker heterozygosity in the population. Ancestral relationships in explain, for instance, genomic relationships in that are not captured by available pedigree; e.g. across nominally unrelated individuals. It will be shown later that this relationship is relative to a population with maximum heterozygosity and is analogous to an fixation index [18].

Christensen [11] estimated the two parameters, and , using maximum likelihood, whereas [1] suggested methods of moments. Closer inspection of Appendix A in [11] leads to the following developments that were not described in Christensen [11] (see “Appendix” of the present paper for more details).

Parameter is such that , with the allele frequency at a random locus . Parameter , with being the number of markers. However, such that if reference alleles are chosen at random across loci, then . From this it follows that:and the genomic relationship matrix is . Interestingly, this matrix is similar to a matrix of IBS relationships, that can be written as:so that (see proof in “Appendix”).

Substituting into the expression gives: , such that for a single population is eight times the variance of allele frequencies in the base population (this variance was described by Cockerham [19]). We stress that to imply that (and ) is a parameter, the variance of allele frequencies across markers [10, 11, 20, 21]. However, can be considered as equivalent to heterozygosity when all markers have an allele frequency of 0.5, that is, the maximum possible heterozygosity.

Multiple populations

In an analogous manner, the relationship between two metafounders and is: , i.e., the covariance across loci between allele frequencies of two populations and . This is almost tautological: the relationship (between two populations in this case) is the covariance between the gene content at a locus. Christensen et al. [6] implicitly show this in Appendix A of their paper. Cockerham [19] and Robertson [22] interpreted 4 as the coancestry between two populations and Fariello et al. [23] used to describe the divergence of populations. Several measures of genetic distance between populations have been developed (e.g. [24]), and most of them contain a term that is related, implicitly or explicitly, to . In particular, the average square of the Euclidean distance can be written as . Thus, .

Estimation

Estimation in a single population

Estimation of is trivial, it is simply half the number of markers. Parameter is proportional to the variance of allele frequencies in the base population. If base population individuals were genotyped, computing allele frequencies and estimating would be trivial. In the next section, we propose methods when this is not the case, i.e. genotyped individuals are related and perhaps several generations away from the base population.

Assuming no pedigree structure i.e. naïve The simplest model assumes that genotyped individuals are unrelated and constitute the base population. For locus , let be a vector of gene contents in the form {0,1,2}, defined as before. The mean of this vector is For each locus, is estimated as the observed mean of , then is computed as the empirical variance across loci of , and because , then and .

Considering pedigree structure At locus , gene content can be seen as a quantitative trait mean of in the base population equal to , where is the allele frequency in the base population and the genetic variance is [25-27]. Cockerham [19] showed that the covariance of gene content of marker between individuals and is a function of their relationship (: . A linear model can therefore be written as:where is an incidence matrix relating individuals in the pedigree to observed genotypes, and is the deviation of each individual from the mean for all individuals [25-27]. Assuming multivariate normality: and .

Equivalently, for the set of genotyped individuals (labelled as “2”), , where is an additive relationship matrix that includes only the genotyped individuals. From this formulation, there are two possible strategies to estimate .

Generalized least squares (GLS) This ignores the prior distribution of and estimates each as a “fixed effect”, using best linear unbiased estimator (BLUE) (or, equivalently, GLS) estimators of separately for each locus. One option is to use the spanning all the pedigree and mixed model equations [25-27]. Equivalently, the corresponding GLS expression is:where is the sum of elements of , and is a weighted sum of genotypes. Then, is estimated as , because , , and it follows that .

Maximum likelihood If allele frequencies in the base population have a distribution, can be considered as drawn from a normal distribution, . Thus is a variance component that can be estimated by maximum likelihood (ML). The equations for given values of and are . An expectation–maximization scheme [28] to obtain ML is as follows. Pick starting values for and . Iterate until convergence on:

For each marker ,

estimate ,

store ,

update as with ;

Update as , where the second part of the expression corresponds to the trace , , the identity matrix, is the relationship structure across levels of and is the prediction error covariance matrix of . As only the diagonal elements of are needed in , its elements can be obtained separately from each single locus analysis.

At convergence, the estimate is . This gives the same estimate as the method based on a Wishart likelihood function [11] with (results not shown).

Estimation in multiple populations

If base populations are considered, the variance component generalizes to , a matrix of variances and covariances between means for marker in population . Across populations, and .

Assuming no pedigree structure

Naïve If relationships across individuals are ignored:where is a matrix, the rows of which sum to 1, and that assigns individuals to fractions of populations, and is a vector with elements for the average of each population. For each locus, can be estimated using least squares and the covariance matrix of across loci gives an estimate of , e.g. for two populations , a two-by-two matrix.

Considering pedigree structure

If there are no crosses between individuals from different populations in the pedigree, the estimation of allele frequencies in each base population can be split in separate analyses by population : , with and the matrix of pedigree-based relationships among individuals in population b, and the analysis proceeds as in a single population. Then, is estimated as the observed matrix of covariances for across loci.

When there are crosses, there are two alternatives.

Generalized least squares (GLS)

The first alternative [27] is to use a genetic groups model [12, 13], as , where contains the fraction of ancestry in individual . This ignores the fact that the variance of gene content, , differs between breeds and crosses. The second, and more exact alternative, is to use the representation where the breeding values are split into within- and across-breed components [29]:with partial relationship matrices for vectors and . The BLUE’s of can be obtained and then estimated as above.

Maximum likelihood (ML)

Analogously to the single population case, an expectation–maximization updated estimate can be obtained using multiple-trait formulations [28], where is the prediction error variance–covariance, e.g. with two populations:

Our implementation of this approach is as follows:

For each marker :

estimate ,

store ,

update as with ;

Update using cross-products within and across populations as, e.g., with two populations:

Step 1 includes an approximation in (1c) because we assume that is the same for all base populations, as in the GLS above, which could be further improved by using partial relationship matrices. This point will be addressed in future research.

Simulation

To assess the quality of genomic predictions using one metafounder, we simulated data using QMSim [30]. The simulation closely followed that in [5] to mimic a dairy cattle selection scheme scenario. A historical population undergoing mutation and drift was generated, followed by a recent population undergoing selection.

First, 100 generations of the historical population were generated with an effective population size of 100 during the first 95 generations, followed by a gradual expansion during the last five generations to an effective population size of 3000. Thirty chromosomes of 100 cM and 40,000 segregating biallelic markers distributed at random along the chromosomes in the first generation of the historical population were simulated. The 40,000 markers were resampled from a larger set of 90,000 markers in order to obtain allelic frequencies from a beta(2,2) distribution, similar to dairy cattle marker data, so that parameter had a true value around 0.40. There were 1500 QTL affecting the phenotype; QTL allele effects were sampled from a Gamma distribution with a shape parameter of 0.4. Mutation rate at the markers (recurrent mutation process) and QTL was assumed to be 2.5 × 10−5 per locus per generation [31]. We used a higher mutation rate than typical (10−8, [32, 33]) to overcome the fact that QMSim is not a coalescent simulator. Phenotypes for a trait recorded only on females with a heritability of 0.30 were simulated.

Then, 10 overlapping generations of selection followed. In each generation, 200 males were mated with 2600 females to produce 2600 offspring by a positive assortative mating design based on EBV. Within the simulation, individuals were selected according to estimated breeding value (EBV) based on pedigree BLUP. In each generation, 40% of males and 20% of females were replaced by selected younger individuals. No restrictions were set to avoid or minimize inbreeding, so highly inbred individuals were found, as a result of strong selection and matings among highly-related individuals. A total of 100 individuals had an inbreeding coefficient higher than 0.20 (mainly found in the last generation), with some individuals having inbreeding coefficients higher than 0.40. True breeding values (TBV) and pedigree information were available for all 10 generations (28,800 individuals in the pedigree), phenotypes were available for all females except in the last generation (14,300 records). The 840 sires of females with phenotypic records were genotyped, as well as 2600 individuals in generation 9 (with records) and 2600 in generation 10 (without records). A total of 20 independent replicates were made. A two-step analysis was carried out using the simulated data. First, we compared several methods to estimate . Then, we tested the quality of genomic predictions using four methods (see section on genomic prediction methods), one of which included one metafounder.

Methods to estimate

Parameter was estimated using four estimation methods. First, the naïve method that does not consider the pedigree structure. Pedigree information was included in three methods: GLS, ML, and the method of moments (MM) in [1]. For a single population, the last method involves estimation of based on summary statistics of (regular pedigree-relationship matrix for genotyped individuals) and (the genomic relationship matrix).

Genomic prediction methods

The EBV of the selection candidates in generation 10 (genotyped and without phenotype records) were estimated using four methods. The first was the pedigree-based BLUP (PBLUP) based on phenotype and pedigree information. The second method was ssGBLUP, in which genomic information is also taken into account. We used the correction of [34] to equate genomic and pedigree average inbreeding and relationships, the default method used in most practical applications [34, 35]. However, the implementation that we used does not include inbreeding in the setup of [36], although it does consider inbreeding in (see below for use of these matrices). The third method was ssGBLUP that includes inbreeding in the setup of and of (ssGBLUP_F). The fourth method was ssGBLUP with the metafounder (ssGBLUP_M), using estimated by GLS since it turned out to be an accurate method to estimate (see the Results section). The four methods used the following inverse relationship matrices: PBLUP: ; ssGBLUP: where is as in [34] and is constructed ignoring inbreeding [36]; ssGBLUP_F: same as ssGBLUP, with correctly constructed; ssGBLUP_M: where with (see the "Methods" section) and is as in [1]. More details are given in “Appendix”. For computation, we used blupf90 [37]. In the case of ssGBLUP_M, we constructed with own software and then used the option user_file in blupf90 (http://nce.ads.uga.edu/wiki/doku.php).

Quality of genomic prediction

Prediction quality was evaluated for all 2600 selection candidates in generation 10. The accuracy of the methods was measured as the Pearson correlation between TBV and EBV. Bias was calculated as the difference between the average TBV and average EBV with respect to the base population (i.e. to the solution of the metafounder for ssGBLUP_M or to 0 for the other methods). Thus, bias is related to estimated genetic progress in the selection candidates. The inflation (often also called bias) of the prediction method was quantified by the coefficient of regression of TBV on EBV. These two statistics correspond to the coefficients and in the Interbull validation method [38], which uses the regression . The mean square error (MSE) of prediction of EBV was calculated as the mean of the squared difference between TBV and EBV. An ideal method should have maximum accuracy, minimum MSE, zero bias, and a regression coefficient of 1. These are not only elegant statistical properties but also have relevance in livestock selection [39-41]. Changes in ranking of the selection candidates were also assessed by calculating the Spearman’s rank correlation coefficients between EBV across methods.

In addition, the quality of variance component estimation was assessed by comparing estimated and simulated heritabilities. For this purpose, variance components were estimated by REML with remlf90 [37] based on the four methods (PBLUP, ssGBLUP, ssGBLUP_F, ssGBLUP_M).

Results

Estimation of

Figure 1 shows boxplots of the differences between the estimates of based on the four methods (MM, Naïve, ML and GLS) and the true values obtained by simulation, for each of the 20 replicates. The simulations were tailored to produce . Methods ML and GLS estimated very accurately. Method MM clearly underestimated , whereas the Naïve method overestimated it. Based on these results, we used estimated by GLS for ssGBLUP_M for prediction. The effect of employing different values of in genomic prediction was assessed to quantify its impact on predictions. Using estimates of based on MM only slightly changed results. For example, the accuracies and slopes of ssGBLUP_M were not affected up to the 4th digit (not shown).

Fig. 1

Differences between estimated and true Gamma, across 20 simulation replicates. Gamma was estimated by generalized least squares (GLS), maximum likelihood (ML), method of moments (MM) and the Naive method

Correlations between TBV and EBV of candidates in generation 10 for each prediction methods are in Table 1 and Fig. 2a. Compared with PBLUP, ssGBLUP_F and ssGBLUP_M increased accuracy by approximately 23 absolute points. This shows an important improvement by including marker information in the prediction and the possibility of generating a small extra gain when also including the metafounder. Method ssGBLUP resulted in a small loss of accuracy compared to ssGBLUP_F and ssGBLUP_M.

Table 1

Accuracy (correlation between TBV and EBV), inflation (regression coefficient of TBV on EBV), bias [average (EBV–TBV)] and mean square error (MSE) for each prediction methods

Prediction method	Accuracy	Inflation	Bias	MSE
PBLUP	0.51 (0.05)	0.98 (0.06)	−0.0003 (0.03)	0.206 (0.01)
ssGBLUP	0.72 (0.03)	0.89 (0.19)	0.2169 (0.04)	0.159 (0.03)
ssGBLUP_F	0.74 (0.02)	0.99 (0.04)	0.1167 (0.04)	0.141 (0.01)
ssGBLUP_M	0.74 (0.02)	0.94 (0.04)	0.0094 (0.03)	0.125 (0.01)

Averages across 20 replicates with standard deviations in parenthesis

Fig. 2

a Correlation of TBV with EVB for each prediction method (accuracy). b Regression slope of TBV on EBV (overdispersion). c Bias [average (EBV–TBV)]

Table 1 and Fig. 2b display the regression coefficient of TBV on EBV, which measures the degree of inflation for each prediction method and should be close to 1. PBLUP and ssGBLUP_F produced values closest to 1. Including genomic data in the prediction based on ssGBLUP resulted in regression coefficients lower than 1, but including the metafounder in ssGBLUP_M gave values closer to 1. Methods ssGBLUP_M and ssGBLUP_F displayed a lower standard deviation compared to the other two methods. Again, method ssGBLUP showed the highest variability.

Biases of EBV obtained with each prediction method are in Table 1 and Fig. 2c. Both PBLUP and ssGBLUP_M were unbiased, whereas ssGBLUP and ssGBLUP_F were biased. The bias was higher for ssGBLUP than for ssGBLUP_F, which was largely due to a single outlier; the median bias was roughly the same for ssGBLUP and ssGBLUP_F. The bias with ssGBLUP_F was equivalent to roughly 0.5 generations of genetic improvement or 0.4 standard genetic deviations. Finally, ssGBLUP_M had the lowest MSE (closer to zero), followed by ssGBLUP_F (Table 1).

Ranking of EBV

The methods were also compared based on rank correlations of EBV with TBV and between methods. A rank correlation of 1 implies that the same candidates would be selected. Results are in Table 2. Rank correlations with TBV were similar to the Pearson correlations in Table 1. Selection decisions differed only slightly when using ssGBLUP, ssGBLUP_F or ssGBLUP_M. Note, however, that this table reports rank correlations among young selection candidates in the last generation and does not address comparisons across generations (e.g. old vs. young animals), which is sensitive to the biases that are reflected in Table 1 [41]. For instance, all young animals would be overestimated by 0.11 with ssGBLUP_F, which results in these young animals looking better than proven sires, which had an accuracy of essentially 1 and no bias. Depending on the selection scheme, this may lead to less than optimal selection decisions.

Table 2

Spearman correlations among TBV and the four EBV for each prediction methods

	EBV PBLUP	EBV ssGBLUP	EBV ssGBLUP_F	EBV ssGBLUP_M
TBV	0.49 (0.06)	0.71 (0.02)	0.72 (0.03)	0.73 (0.02)
EBV PBLUP		0.56 (0.05)	0.62 (0.04)	0.64 (0.04)
EBV ssGBLUP			0.99 (0.01)	0.98 (0.01)
EBV ssGBLUP_F				0.99 (0.002)

Averages across 20 replicates with standard deviations in parenthesis

Estimation of variance components

Figure 3 shows estimates of heritability obtained with three of the four methods (PBLUP, ssGBLUP_F and ssGBLUP_M). The estimates obtained using ssGBLUP did not converge for six of the 20 simulation replicates. Convergence was achieved in those cases by weighting the submatrix in by instead of 1 [42] but poor quality estimates were obtained and are, therefore, not reported.

Fig. 3

Estimated heritability for PBLUP, SSGBLUP_F and SSGBLUP_M considering the 20 replicates. The dotted line shows the simulated heritability of 0.30

Estimates were generally lower than the simulated true heritability (0.30). The lowest estimates were obtained with ssGBLUP_F. Including the metafounder improved estimates compared to ssGBLUP_F and reduced variability of estimates compared to PBLUP.

Discussion

In this work, we have addressed the complex issue of conciliation of marker and pedigree information in genetic evaluation. Powell et al. [43] argued that both IBS (at the markers) and identity-by-descent (IBD) are compatible notions because they are both measures of identity at causal genes. However, incompatibility appears when mixing both types of relationships [5, 34, 44, 45]. Legarra [7] suggested that, in order to compare genetic variance across IBD, IBS or other measures of relationships, a common reference must be chosen. Similar (but not identical) to [43], in this work we used a fixed reference ( constructed as a cross-product of {−1,0,1} genotypic codes) and tailored (IBD, pedigree) to fit (IBS, markers). Compared to previous approaches, using a fixed reference has the advantage that genomic relationships are immutable (i.e. adding more genotyped individuals to the database does not change the existing relationships) and they do not depend on pedigree depth, which by construction is always limited and, in animal breeding, often heterogeneous. Our approach is in fact very similar to using IBS as measure of identity. We used a genomic relationship matrix , whereas the matrix of IBS is (see proof in “Appendix”). In GBLUP with associated variance component estimation, when all animals are genotyped, using a model with or the matrix proposed here yields identical EBV, as the term ½ in gets absorbed into the variance component and the constant gets absorbed into the fixed part of the linear mixed model [7, 46]. However, matrix rather than must be used in ssGBLUP_M because is not compatible with pedigree relationships. In [4], the (fixed effect) intercept term models, identical to [5], the difference between genetic values of individuals in the base and genotyped individuals. These intercept terms play therefore a similar role as metafounders.

Easy estimation of ancestral relationships

Derivations in the Theory section show that estimation of ancestral relationships based on (one base population) and (several base populations) can be framed within the classic linear model approach of quantitative genetics [19], which has recently been used for gene content [14, 25–27]. This approach is easy to understand and compute. Also, can be understood, just like heritability, as an unobserved base population parameter that does not change with additional data (although its estimate may change). Therefore, an accurate estimate of can be used repeatedly without the need for re-estimation, as is customary in livestock genetic evaluation. This contrasts with “centering” of marker covariates, which changes with every new genotype. If all base allele frequencies were known exactly, then there should be no need to use metafounders, as relationship matrices can be appropriately constructed [14].

In the work presented here, the simplest methods (Naïve and method of moments) yielded biased (upwards and downwards, respectively) estimates of ; the naïve method because it ignores that allele frequencies tend to drift to their extreme values as generations progress, and the method of moments because it implicitly assumes that genotyped individuals are a random sample from a particular generation.

Equivalence of ancestral relationships with second moments of allele frequencies also shows a strong relation with population genetics theory, which will be detailed in the next paragraph.

Relationship between metafounders parameter and fixation index

The fixation index [17] is a measure of diversity of a set of populations with respect to a reference population, usually the pool of all populations. In this approach, each population is assumed to be a random sample from all possible populations that could be sampled according to the evolutionary process described by . Conceptually, is a parameter to be estimated [18, 19], and it is not a statistic computed from the data. The usual definition of for a particular biallelic locus is:where is the variance of allelic frequencies across populations and is the allele frequency of the conceptual combined population. If we consider that the variance of allele frequencies applies across loci and not across populations, it follows that because reference alleles are taken at random. In this case:

Our interpretation of this link between and is as follows. Jacquard [47] called the “inbreeding coefficient of a population”. Cockerham [16] modelled as an intraclass correlation, “the coancestry of the line with itself”, in other words, the probability that two gametes taken at random from the population are identical. Thus, it makes perfect sense to consider that the additive relationship (which is twice the coancestry value) of a group with itself is . This is the interpretation of the coefficient in Legarra et al. [1]. Note that the assumption that is automatically fulfilled if reference alleles are chosen at random across loci (i.e., they are neither the most frequent nor the least observed allele).

Alternatively, [1] showed that for a population with self-relationships equal to , the average heterozygosity is , i.e. the variance is reduced by an amount equal to from the conceptual population with heterozygosity 1. Thus can be interpreted as if the is taken as a measure of homozygosity.

Consequences of using metafounders in genomic evaluation

Genomic estimates of breeding values are invariant to allele coding [46] when all individuals are genotyped. However, this is not the case when pedigree and marker informations are combined, as in ssGBLUP. In this work, we have shown that, even in the presence of complete pedigree and a single base population, use of metafounders in ssGBLUP_M leads to slightly more inflated and less biased EBV, lower MSE, and nearly unbiased estimates of heritability compared to ssGBLUP_F. Bias, defined as E(EBV–TBV), is typically overlooked in genomic prediction, but in an example of biased evaluation, Henderson [48] recognized that “sires of later generations appeared to be under-evaluated relative to older sires”. Overdispersion, also called bias in recent literature (e.g. [38]), may also have a dramatic impact in practice [39-41] and the trade-off between bias and variance needs further study. For instance, Vitezica et al. [5] found that ssGBLUP_F was unbiased but had some overdispersion, which likely depends on the data structure, including which individuals are genotyped.

In addition, use of metafounders allows a clear definition of genomic relationships because relationships do not depend on pedigree depth or completeness or on changes in allele frequencies as new data is added. In addition, a high-dimensional parameter (i.e. base allele frequencies) is substituted by a low-dimensional one (matrix ).

The poor performance of ssGBLUP compared to ssGBLUP_F is likely due to the presence of highly inbred individuals because ssGBLUP ignored inbreeding in the setup of . This relates to the interpretation of parameter , as used in early studies of ssGBLUP [42]. An application of ssGBLUP for type traits in Holstein [42] experienced convergence problems, which were eliminated when was multiplied by and which increased accuracy of predictions. However, the nature of parameter was not known [49]. In those studies, the inverse of the numerator relationship matrix was constructed using Henderson’s rules while ignoring inbreeding [36], while the submatrix included inbreeding. As a result, the elements in the latter matrix were too large. In addition, genotyped animals were on average unrelated in but not in , which can be corrected by scaling , as in [5]. However, this resulted in the elements in to be too large for younger animals relative to . Both these problems are partially circumvented by putting a weight on . When was constructed while considering inbreeding, the optimal value of in an analysis of Holstein dairy cattle increased from 0.7 to 0.9 (Masuda, personal communication, 2016). However, the metafounder approach provides a more principled solution to this problem. Also, following these experiences, should always be constructed while considering inbreeding to avoid infrequent but pathological problems.

Conclusions

Metafounders have relationships that are closely related to fixation indices and proportional to covariances of allele frequencies in base populations. Use of metafounders can be simplified by new methods to estimate the covariance of base allele frequencies. We verified by simulation of a selected population that, in a single population, both GLS and ML are unbiased and computationally efficient. In the same simulation, use of metafounders in ssGBLUP led to more accurate and less biased evaluations, and also to more accurate estimates of genetic parameters. We propose a genomic relationship matrix that refers to a population with ideal base allele frequencies equal to 0.5. This matrix is similar to an IBS relationship matrix (up to scale factors), does not change with new data, and is compatible with pedigree data if metafounders are used. In the simulated data, pedigrees were perfectly known. Future work with real datasets in more complex settings—purebreds and their crosses [50, 51], and selected populations with unknown parent groups [13] will investigate the feasibility and accuracy, in practice, of using metafounders in ssGBLUP.

	AA	Aa	aa
AA	2	1	0
Aa	1	1	1
aa	0	1	2

	AA	Aa	aa
AA	1	0	−1
Aa	0	0	0
Aa	−1	0	1

No parent known	1
One parent known	0.75-F_known4^-1
Two parents known	0.5-F_sire4-F_dam4^-1