A comparison of marker-based estimators of inbreeding and inbreeding depression

Background The availability of genome-wide marker data allows estimation of inbreeding coefficients (F, the probability of identity-by-descent, IBD) and, in turn, estimation of the rate of inbreeding depression (ΔID). We investigated, by computer simulations, the accuracy of the most popular estimators of inbreeding based on molecular markers when computing F and ΔID in populations under random mating, equalization of parental contributions, and artificially selected populations. We assessed estimators described by Li and Horvitz (FLH1 and FLH2), VanRaden (FVR1 and FVR2), Yang and colleagues (FYA1 and FYA2), marker homozygosity (FHOM), runs of homozygosity (FROH) and estimates based on pedigree (FPED) in comparison with estimates obtained from IBD measures (FIBD). Results If the allele frequencies of a base population taken as a reference for the computation of inbreeding are known, all estimators based on marker allele frequencies are highly correlated with FIBD and provide accurate estimates of the mean ΔID. If base population allele frequencies are unknown and current frequencies are used in the estimations, the largest correlation with FIBD is generally obtained by FLH1 and the best estimator of ΔID is FYA2. The estimators FVR2 and FLH2 have the poorest performance in most scenarios. The assumption that base population allele frequencies are equal to 0.5 results in very biased estimates of the average inbreeding coefficient but they are highly correlated with FIBD and give relatively good estimates of ΔID. Estimates obtained directly from marker homozygosity (FHOM) substantially overestimated ΔID. Estimates based on runs of homozygosity (FROH) provide accurate estimates of inbreeding and ΔID. Finally, estimates based on pedigree (FPED) show a lower correlation with FIBD than molecular estimators but provide rather accurate estimates of ΔID. An analysis of data from a pig population supports the main findings of the simulations. Conclusions When base population allele frequencies are known, all marker-allele frequency-based estimators of inbreeding coefficients generally show a high correlation with FIBD and provide good estimates of ΔID. When base population allele frequencies are unknown, FLH1 is the marker frequency-based estimator that is most correlated with FIBD, and FYA2 provides the most accurate estimates of ΔID. Estimates from FROH are also very precise in most scenarios. The estimators FVR2 and FLH2 have the poorest performances. Supplementary Information The online version contains supplementary material available at 10.1186/s12711-022-00772-0.


Results:
If the allele frequencies of a base population taken as a reference for the computation of inbreeding are known, all estimators based on marker allele frequencies are highly correlated with F IBD and provide accurate estimates of the mean ΔID. If base population allele frequencies are unknown and current frequencies are used in the estimations, the largest correlation with F IBD is generally obtained by F LH1 and the best estimator of ΔID is F YA2 . The estimators F VR2 and F LH2 have the poorest performance in most scenarios. The assumption that base population allele frequencies are equal to 0.5 results in very biased estimates of the average inbreeding coefficient but they are highly correlated with F IBD and give relatively good estimates of ΔID. Estimates obtained directly from marker homozygosity (F HOM ) substantially overestimated ΔID. Estimates based on runs of homozygosity (F ROH ) provide accurate estimates of inbreeding and ΔID. Finally, estimates based on pedigree (F PED ) show a lower correlation with F IBD than molecular estimators but provide rather accurate estimates of ΔID. An analysis of data from a pig population supports the main findings of the simulations.
Conclusions: When base population allele frequencies are known, all marker-allele frequency-based estimators of inbreeding coefficients generally show a high correlation with F IBD and provide good estimates of ΔID. When base population allele frequencies are unknown, F LH1 is the marker frequency-based estimator that is most correlated with F IBD , and F YA2 provides the most accurate estimates of ΔID. Estimates from F ROH are also very precise in most scenarios. The estimators F VR2 and F LH2 have the poorest performances.

Background
The characterization and management of inbreeding are important topics in conservation [1] and in population and quantitative genetics [2][3][4]. Proper estimation of individual inbreeding coefficients is necessary to estimate inbreeding depression, i.e., the change in the mean of a quantitative trait, which mainly results from the expression of recessive alleles with their effect being hidden in heterozygotes and exposed in homozygotes by inbreeding (see [2], Chapter 10). The inbreeding coefficient of an individual is defined as the correlation between homologous alleles [5] or the probability of identity-bydescent (IBD) of homologous alleles [6], and has been traditionally obtained from pedigree data [7]. However, inbreeding can also be obtained from genomic information (see e.g., [8]) and the development of high-density single nucleotide polymorphism (SNP) panels for an increasing number of species gives the opportunity to develop and implement methods that can accurately estimate the true realized inbreeding [9]. Genomic information also enables the levels of genetic variation across the genome to be evaluated [10][11][12] or the identification of genomic regions that are responsible for inbreeding depression [13,14].
There are many metrics to estimate relatedness and inbreeding from marker data that are applied on a SNPby-SNP basis and need the knowledge of marker allele frequencies [15][16][17][18][19][20][21][22][23][24]. Another useful marker-based measure to estimate inbreeding is the proportion of the autosomal genome that is composed of runs of homozygosity (ROH) [25], as well as other alternative analogous metrics [26]. ROH measures of inbreeding provide estimates of IBD and have been shown to be useful estimators of the inbreeding coefficient and as a basis for estimating inbreeding depression [9,27,28] but their performance depends on the criteria used to define a ROH [29,30]. Marker frequency-based measures of inbreeding are expected to provide unbiased estimates of the average IBD measures of relatedness and inbreeding relative to an initial generation or base population that is assumed to consist of non-inbred and unrelated individuals, when the allele frequencies in that initial generation are known [19,22,31,32]. However, allele frequencies in the base population are usually not available and inferring them is difficult [33], thus frequencies in the current generation are often used instead. In that case, marker frequencybased measures of inbreeding provide deviations from Hardy-Weinberg proportions or correlations of allele frequencies of homologous genes relative to the current generation of the population rather than estimates of IBD [22], and the estimates are confounded by the average coancestry of individuals with other individuals in the sample taken for analysis [34].
The different measures of molecular inbreeding have been evaluated in multiple studies, both based on empirical data [13,14,24,[34][35][36][37][38][39][40][41][42][43][44], and based on simulated data or theoretical analyses [9,32,34,[45][46][47]. Some empirical studies have used inbreeding estimators based on known [12] or estimated [35] allele frequencies from a predefined base population. In other studies, a frequency of 0.5 is assumed for all SNP alleles, which results in inbreeding estimates that have intermediate to high correlations with pedigree-based estimates or IBD [35,37,48]. Nevertheless, in most studies, allele frequencies in the current population are used to obtain different marker-based inbreeding estimates and these are compared with each other and with pedigree-based inbreeding estimates, with rather contrasting results. For example, the correlation of coefficients of inbreeding obtained from pedigrees with those obtained from markers is often intermediate ranging from 0.4 to 0.8 [13,14,39,41,42,49], but can sometimes be low or even negative [39,41,44,50,51].
Previous simulation studies have investigated the accuracy of different molecular measures of inbreeding to estimate IBD. One of the most comprehensive studies is that from Wang [23], who showed that when the number of markers is large enough (around 10,000), markerbased estimates of inbreeding are more correlated with IBD than pedigree-based estimates. A later simulation study by Forutan et al. [48], who extended the analyses to artificially-selected populations, showed that markerbased inbreeding coefficients were close to the IBD values when the frequencies of the base population are known, as expected. However, when constant allele frequencies of 0.5 were used, marker-based coefficients overestimated IBD, although the correlation of markerbased estimates with IBD values was close to 1.
The accuracy of using molecular measures of inbreeding based on current allele frequencies to estimate the rate of inbreeding depression (ΔID) has also been addressed by different studies [9,23,32,46,47]. These studies have shown that marker-based measures of inbreeding and relatedness are more precise and more powerful to detect inbreeding depression than pedigreebased measures. In addition, simulation results have indicated that some estimators of molecular inbreeding, such as those obtained from the correlation of uniting gametes [20] and those from ROH, generally provide good estimates of ΔID but can give biased estimates in some cases. However, comprehensive studies are still needed to elucidate the performance of the use of different markerbased estimators to estimate ΔID under some scenarios with non-random contributions from parents to progeny. For example, in conservation programs it is common to equalize contributions from parents to offspring in order to maintain the highest genetic diversity [1]. In addition, many populations are subjected to artificial selection, where contributions from parents to progeny can be far from random.
The objective of this study was to evaluate how different estimators of molecular inbreeding perform as measures of IBD and how precise they are when estimating ΔID, with or without knowledge of allele frequencies in the base population of reference. We focus on the most popular marker-based estimators that make use of SNP allele frequencies or ROH [42][43][44][51][52][53][54][55][56][57][58] and that are readily obtained by commonly used software such as PLINK [59,60] and GCTA [36]. We carried out computer simulations that assumed a model of selection against deleterious mutations, with random or equal contributions from parents to progeny, or artificial selection for a neutral quantitative trait. We also analysed a pedigreed population of Iberian pigs where equalization of family sizes has been intended over 23 generations, and compared values of inbreeding and inbreeding depression obtained for the different measures of inbreeding based on a simulation of this population with their observed empirical values.

Simulation of populations under different breeding program designs Large initial population
We used the software SLiM 3 [61] to simulate the genome of a large population of a sexually reproducing diploid species with random mating and discrete generations. As a model, we assumed the genomic characteristics of the pig genome [62], as some of the simulation results are compared with those obtained from empirical data of a strain of Iberian pigs (the Guadyerbas strain). A population with a constant size of N = 500 individuals was simulated for 5000 generations, which is the approximate effective population size estimated for the ancestral population of the Guadyerbas strain [63]. We assumed a total genome size of 2250 Mb, with 18 chromosomes of 125 Mb each. A summary of the main parameters used in the simulations is in Table 1. A recombination rate of 8 × 10 -9 was assumed, with crossover frequencies following a Poisson distribution with a uniform probability across the genome and without interference. This implies a genome length of about 1 Morgan per chromosome and a recombination rate of 0.8 cM/Mb, which is close to that found in the pig genome [62]. The mutation rate per nucleotide was assumed to be 4 × 10 -9 in order to obtain a number of SNPs close to that available for the pig population, but a model with a higher mutation rate (9 × 10 -9 ) was also considered to investigate a higher density of SNPs.
In order to evaluate ΔID for fitness, it was assumed that 5% of the mutations were deleterious. The fitness of the wild type genotype, the heterozygote, and the mutant homozygote were 1, 1 − sh , and 1 − s , respectively, where the homozygous selection coefficient s was obtained from an exponential distribution with mean s = 0.1, and the dominance coefficient h was obtained from a uniform distribution between 0 and e (−ks) , where k was set to obtain an average value of h = 0.2 [4,64]. Individual genotypic fitness values were obtained multiplicatively across loci. We also assumed an alternative model with a lower mean effect of deleterious mutations ( s = 0.025) and a higher mutation rate ( Table 1). The distribution of mutation effects and dominance coefficients assumed in the two models are shown in Additional file 1: Fig. S1. Table 1 Simulation parameters for the default model, a model with a higher density of SNPs, and an alternative model of deleterious mutations with a lower mean effect than for the default model

Small populations with random contributions, equalization of parents´ contributions, or artificial selection
An in-house C program was used to simulate small populations of size N = 20 or 100 breeding individuals over 10 (for N = 20) or 50 (for N = 100) discrete generations, such that the average inbreeding coefficient ( F ) in the last generation was about the same for both population sizes. Populations of size N = 20 and N = 100 were also simulated for 20 and 100 generations, respectively, in order to approximately double the average inbreeding coefficient in the last generation. To start these populations, individuals were randomly sampled from the large initial population of 500 individuals that was simulated as described previously. The first generation ( t = 0) was considered to be the base population that all inbreeding coefficients refer to. In one set of simulations, generation t = 10 was considered the base population. The genome length, the number of chromosomes, the recombination rate between nucleotides, and the fitness effects were as described previously. It was assumed that the same SNP panel was used to compute measures of inbreeding in the initial and final generations. Thus, only SNPs that segregated at generation t = 0 were considered in the analyses and the mutation rate was set to zero for the subsequent generations. In order to estimate IBD probabilities, about 4500 multiallelic loci uniformly distributed across the genome were also simulated, with unique alleles assigned to each individual in generation t = 0. For the random contributions (RC) and equalization of parents´ contributions (EC) scenarios, a number N of individuals were randomly taken from the large initial population to create each replicated line. For the RC scheme, parents were randomly mated to generate N progeny, allowing polygamy but avoiding self-fertilization. For the EC scheme, the same procedure was followed, except that, in each generation, the N individuals were randomly paired in couples and two offspring were obtained from each couple. For the artificial selection scheme (SEL), T individuals (100 or 500) were taken at random from the large initial population to generate the base population of each replicate. From these individuals, N (20 or 100, respectively) were artificially selected based on their own phenotype for a quantitative trait (described below), implying a selected proportion of 20%. Two thousand replicates were carried out for each simulated scenario.
The genetic basis of the quantitative trait was established as described in the following. From the SNPs of the large initial population, a random 5% were chosen as quantitative trait loci (QTL) and an effect on the quantitative trait was assigned to one of the two alleles. Effects were taken from a normal distribution with mean zero and standard deviation 0.1, assuming additive gene action. Additivity was also assumed between loci, and thus genotypic values were obtained as the sum of the effects of all QTL. The phenotypic value was obtained by adding an environmental deviation to the genotypic value, which was taken from a normal distribution with mean 0 and environmental variance 2, a value chosen such that it generated an initial heritability for the trait of around 0.5. The quantitative trait was assumed to be not related to fitness.
Pedigrees were recorded over generations to obtain the pedigree inbreeding coefficient of each individual in the last generation ( t = 10, 20, 50 or 100) to obtain an estimate of F PED . The identity-by-descent inbreeding coefficient ( F IBD ) in the last generation was computed using the multiallelic loci. The genotypes of the SNPs (about 25,000-40,000 in the default model and more than twice as many in the higher density model; Table 1) were used to obtain the molecular inbreeding coefficients for each individual in the last generation. Loci that affected the selected quantitative trait or fitness were excluded when computing any measure of inbreeding. Marker allelefrequency based measures of F (described below) were obtained using the known allele frequencies of the initial generation ( t = 0, considered the base population of reference), those of the last (current) generation ( t = 10 or 20 for N = 20 and t = 50 or 100 for N = 100), or set to a constant value of 0.5. Runs of homozygosity were also obtained for individuals in the last generation. The mean and variance of the F values of individuals in the last generation were obtained for each estimator of inbreeding. Pairwise Pearson's correlations were also obtained between all inbreeding measures. Average estimates of inbreeding from the different molecular estimators were compared with the inbreeding coefficient obtained based on IBD values ( F IBD ) in the last generation.
Estimates of the rate of inbreeding depression (ΔID) were obtained for fitness as the regression coefficient of the logarithm of individual fitnesses in the last generation on the corresponding estimated inbreeding coefficients (see [2], p. 262). Additional file 2: Table S1 and Additional file 3: Fig. S2 illustrate the removal of deleterious mutations by genetic purging and drift, the reduction in the expected inbreeding depression in consecutive generations, and the change in additive and dominance variance for fitness for one of the simulated scenarios.
Power to detect inbreeding depression was obtained following Wang [23], by counting the percentage of replicates for which the value of ΔID was significantly different from zero with a probability P < 0.05. A bootstrap method based on 1000 resampling sets with replacement of the 2000 pairs of simulated estimates was used to compare the estimates of ΔID obtained from F IBD with those from the different molecular marker estimators. To confirm the bootstrap results, a paired sample t-test was also carried out using R Commander [65].

Simulation of the Iberian pig population
The Guadyerbas population consists of a pedigreed population with 1206 individuals born across 23 cohorts, from which 219 (those born in cohorts 17-23) were genotyped with the Illumina Porcine SNP60 Bead Chip, which contains 35,519 SNPs that segregate in the Iberian breed [66]. Phenotypes for litter size were available for 832 females (2712 litters, with a mean litter size of 7.4 piglets and standard deviation of 2.3) that were born across all cohorts (i.e., from 0 to 23), so that ΔID could be calculated by including F PED as a covariate in an animal model. The detailed statistical model used is given by Saura et al. [14]. Estimates of ΔID from molecular measures of F were also obtained for 103 females (250 litters, with mean litter size of 7.1 piglets and standard deviation of 2.4) that were genotyped in cohorts 17-23 [14].
We carried out a gene dropping simulation of the whole known pedigree and ascribed to the individuals of the initial cohort the simulated genotypic data from the large initial population of 500 individuals, as described before. The pedigree was run 2000 times, starting from random samples from the large initial population to generate the genotypic values for the individuals of the initial cohort (cohort 0). The length of the genome, the number of chromosomes and the recombination rate assumed between nucleotides were the same as those detailed in the previous sections. Simulated offspring were obtained for each mating following the pedigree in the absence of selection and mutation. The different measures of F were obtained for all individuals across cohorts. The mean and variance of the different estimators of inbreeding and the correlations between them were calculated for the last cohort (23), either assuming the allele frequencies of cohort 17 (considered the base population to which F estimates refer to), assuming the frequencies of the current generation (cohort 23), or assuming frequencies equal to 0.5. Finally, using the available data for litter size from 103 females (born in cohorts 17-23), ΔID was obtained for each molecular F estimator using the simulated genotypes. The simulation results were compared with the empirical estimates.

Molecular measures of the inbreeding coefficient
The following estimators of the inbreeding coefficient were calculated based on the SNP data: where S is the total number of markers, x k is the number of minor alleles of marker k (i.e., 0, 1 or 2 copies), and p k is the frequency (initial, current or 0.5) of the minor allele. This estimator was proposed by VanRaden [19] and it is based on the variance of additive genetic values. a metric proposed by Amin et al. [67], as cited by Van-Raden [19]. It only differs from F VR1 in that the weighting by the variance of allele frequencies is done for each SNP rather than by the summation of variances for all SNPs, so that it gives a larger weight to rare alleles.
which is a coefficient based on the correlation between uniting gametes [20], such that homozygous genotypes are weighted by the inverse of their allele frequencies, and it is expected to show lower sampling variance than other estimators [36].
an analogous estimator to that of Yang et al. [20] ( F YA2 ) but considering the summation of terms separately in the numerator and denominator, as recommended by Zhang et al. [34]. the metric that was initially proposed by Li and Horvitz [16], which gives the deviation of the observed frequency of homozygotes from the expected values under Hardy-Weinberg equilibrium when current allele frequencies are considered.
for which the different weighting in relation to F LH1 is analogous to that between F VR2 and F VR1 , and between F YA2 and F YA1 .
which is obtained by setting p k = 0.5 for any of the previous estimators.
which is the average number of homozygous SNPs, which has a correlation of 1 with F LH 1 and F q05 . The above estimators may receive different names in different studies and a summary of this nomenclature is in Table 1 of Villanueva et al. [12]. Here, we used a terminology, as simple as possible, by identifying the abbreviation for the authors who originally proposed the estimators (VR, YA, LH) plus a suffix 1 to indicate that the measure is obtained from a ratio of averages over loci (VR1, YA1, LH1) and a suffix 2 when the same estimator is obtained from the average of the ratios (VR2, YA2, LH2). This may simplify the existing "Babel tower" [44] of nomenclature.
The estimators are based either on genetic drift derivations ( F VR1−2 and F YA1−2 ) or on homozygosity derivations ( F LH 1−2 and F HOM ) and the equivalence between both approaches was already demonstrated by Cockerham [68]. A summary of the relationships between the estimators is given in the Supplementary Appendix of Caballero et al. [32], and a derivation of F YA2 from the correlation between uniting gametes is given in the Appendix of this paper.
The measures F VR2 , F YA2 , F LH 1 and F LH 2 can readily be obtained by the software PLINK and GCTA ( F I , F III , F and F II , respectively), while F VR1 can be obtained from the diagonal of the genomic relationship matrix obtained by GCTA with the option-make-grm-alg 1.
Note that, for these molecular measures to provide estimates of F IBD with reference to a given base population, the calculations must be done using the allele frequencies ( p k ) of all SNPs that segregate in the base population, rather than of only SNPs that segregate in the generation for which the estimate of inbreeding is computed. Estimates F HOM and the coefficient obtained assuming frequencies equal to 0.5 ( F q05 ) were obtained using only SNPs that are segregating in the current generation.
Finally, estimates of the inbreeding coefficient based on ROH ( F ROH ) were obtained as: where L ROH is the sum of the lengths of all ROH across the genome of the individual, and L is the genome length [25]. We used the software PLINK [59,60] with default options, i.e., a minimum of 100 SNPs per ROH (30 SNPs in the simulated pig population to match the conditions assumed by Saura et al. [14]), at least 1 SNP per 50 kb (100 kb in the simulated pig population), a scanning window of 50 SNPs, and two SNPs in the same ROH no more than 1000 kb apart. We considered ROH lengths longer than 1 Mb ( F ROH −1 ) or 5 Mb ( F ROH −5 ).

Comparison of estimators of inbreeding and inbreeding depression
Default model The average F IBD under EC was about one half of that under random contributions, and that under artificial selection was slightly larger than that without selection, as expected ( Fig. 1, first row). The average F IBD was well estimated by F PED , and also by the five markerbased inbreeding coefficients when allele frequencies of the base population are known (blue bars), with a slight underestimation due to the excess of heterozygotes expected for markers in finite populations . When allele frequencies in the current generation were used instead (yellow bars), all estimators gave the expected average deviation from Hardy-Weinberg proportions (values slightly below zero, expected to be around − 1/2N). Estimates from direct average homozygosity of SNPs ( F HOM ) gave large overestimations of F IBD , as expected. The average values of F ROH were slight overestimations of F IBD , particularly when including shorter fragments ( F ROH −1 ). Using allele frequencies equal to 0.5 when computing the marker-based estimators ( F q05 , purple bars) gave also substantial overestimations of F IBD .
The variances between individual estimates of inbreeding were substantially smaller for F PED , F HOM and F q05 than for F IBD (Fig. 1, second row, VF). When the initial allele frequencies were known (blue bars), variances of molecular estimators of inbreeding were very similar to the variance of F IBD , except that for F VR2 (particularly for the SEL scheme), which had the largest variance. When allele frequencies of the current generation were used (yellow bars), the variances for molecular estimators of inbreeding increased substantially, although not so much in the EC scenario. The estimator that had the closest variance to that of F IBD was F YA2 . The variances of F ROH were also close to those of F IBD .    The third row of the graphs in Fig. 1 shows the correlation between the estimated and true F IBD values. When the initial allele frequencies were known (blue bars), all molecular estimators had a very high correlation with F IBD , except for F VR2 . The F ROH estimators also had a very high correlation with F IBD , while the F PED estimator had the lowest correlation with F IBD . When using allele frequencies in the current generation (yellow bars), F LH 1 , F HOM , and F q05 (which have a correlation of 1 among them) generally showed the highest correlations with F IBD .
The fourth row of the graphs in Fig. 1 shows estimates of the inbreeding depression rate (ΔID) compared to the estimate based on F IBD (horizontal dotted line). When base population frequencies were known (blue bars), all measures of inbreeding provided estimates of inbreeding depression that were close to those obtained with F IBD , with estimates based on F VR2 and F LH 2 slightly biased downwards, particularly in the SEL scenario. Estimates of inbreeding depression based on F PED and F q05 were also rather accurate but F q05 provided a substantial overestimation of inbreeding depression under the EC scenario. When frequencies in the current generation were used (yellow bars), F YA2 gave the most accurate estimates of inbreeding depression for the three scenarios, while F HOM resulted in very large overestimations and F VR2 and F LH 2 in the largest underestimations. Estimates of ΔID based on F ROH were very close to those from F IBD . Statistical tests to compare the mean difference between the estimates of ΔID based on F IBD with those based on the different marker-based estimators are shown in Additional file 4: Table S2. All differences were significantly greater than 0 except for estimates based on F ROH −1 in the RC and EC scenarios, based on F PED in the EC scenario, and based on F YA2 when using the current frequencies in the EC scenario.
Additional file 5: Fig. S3 shows results analogous to those in Fig. 1 but after 20 instead of 10 generations, so that the average inbreeding coefficient was doubled, while Additional file 6: Fig. S4 shows results for N = 100 run after 50 generations. In general, the results in these two figures were similar to those presented in Fig. 1, although there were some relevant differences. With N = 20 after 20 generations [see Additional file 5: Fig. S3], using base population allele frequencies equal to 0.5 (purple bars) resulted in underestimation of average inbreeding and ΔID for the RC and SEL schemes. With N = 100 after 50 generations [see Additional file 6: Fig. S4], the increase of F by selection (scheme SEL) was more evident than in the previous cases (with some overestimation of the molecular estimates), and the estimates of F PED were more clearly biased downwards. Additional file 7: Table S3 gives the minimum and maximum values of individual inbreeding coefficients for the simulations presented in Fig. 1 and in Additional file 5: Fig. S3 and Additional file 6: Fig. S4. As expected, F IBD , F PED , F HOM and F ROH , which are probabilities, ranged from 0 to 1, while F VR1−2 , F YA1−2 and F LH 1−2 , which are correlations or deviations from Hardy-Weinberg proportions, ranged from -1 to 1. However, some values of F VR2 were greater than 1 and some values of F LH 2 were smaller than -1.
The power of the different measures of inbreeding to detect inbreeding depression is shown in Additional file 8: Fig. S5 for the scenarios considered in Fig. 1 and in Additional file 5: Fig. S3 and Additional file 6: Fig. S4. When allele frequencies of the base population were known, all estimators (as well as F ROH and F q05 ) resulted in approximately the same power as F IBD , although some differences were observed between estimators in the case with N = 100 [see Additional file 8: Fig. S5c]. When allele frequencies in the base population were unknown, the most powerful estimator was generally F YA2 (or F YA1 , with a slightly higher power in some occasions). The F VR2 and F LH 2 estimators generally showed the lowest power. The power of F PED to detect inbreeding depression was usually lower than that of marker-based estimators.
Correlations between individual values for all estimators of inbreeding are shown in Additional file 9: Fig.  S6. When allele frequencies in the base population were known [see Additional file 9: Fig. S6a-c], pairwise correlations between marker-based estimators were generally very large (≥ 0.9), with the correlation between F VR2 and F LH 1−2 being somewhat lower. Correlations of F PED with the marker-based estimators were lower than correlations among the latter estimators. When the allele frequencies in the current generation were used, correlations were generally lower than when frequencies in the base population were used but still moderately high [see Additional file 9: Fig. S6d-f ]. Some correlations with F VR2 or F LH 2 were close to zero or negative. The F ROH estimator was highly correlated with the other estimators except with F VR2 , particularly in the SEL scenario.

Alternative models
In order to assess the robustness of the previous results, we also ran scenario RC under alternative conditions. Additional file 10: Fig. S7a shows results for the scenario with N = 20 after 10 generations but with double the density of SNPs than in Fig. 1 (see Table 1). Additional file 10: Fig. S7b shows results for the same scenario except that the effect of deleterious mutations was 1/4 of that used for the default model. In addition, Additional file 10: Fig. S7c presents results for N = 100, as in Additional file 6: Fig. S4, but after 100 generations instead of 50. In general, results for these alternative models were very similar to those of the default model, confirming the main results. Finally, Additional file 10: Fig. S7d shows that results were similar to those of Additional file 5: Fig. S3 with N = 20 after 20 generations but with the base population set up at generation 10 rather than generation 0. The only relevant difference was that F ROH overestimated F IBD and ΔID, which is expected because F IBD is relative to generation 10 whereas F ROH provides values of inbreeding relative to more ancestral generations.
Because the general patterns observed for the different scenarios (RC, EC, SEL, and different N values) were rather similar, Fig. 2 (Fig. 2a). When the allele frequencies of the base population were known (blue bars), all molecular estimators of inbreeding showed a good accuracy for estimating F IBD , although F VR2 was slightly less accurate. When allele frequencies in the current generation were used (yellow bars), F LH 1 (or F HOM ) gave the most accurate estimation of F IBD . Good estimates of F IBD were also obtained when the frequencies of the base population were assumed to be 0.5 (purple bar) or when F ROH estimates were considered (brown bars). Figure 2b gives the average difference between the estimate of ΔID based on the different estimators compared to the estimate based on F IBD for all scenarios. Positive values imply underestimation of ΔID and negative values imply overestimation. Both F PED and all molecular estimators, except F VR2 and F LH 2 , gave relatively good average estimates of ΔID when allele frequencies in the base population were known or when a frequency of 0.5 was assumed. If allele frequencies in the base population were not known, the best estimator of ΔID was F YA2 , although F ROH also provided little biased estimates of ΔID. The F HOM estimator showed large overestimations of ΔID. Estimates were obtained assuming allele frequencies from cohort 17, the current generation (cohort 23), or equal frequencies of 0.5. The figure also shows simulated values of ΔID (Fig. 3e, f ). The simulation results generally agreed with those obtained for the discrete generation simulations, although results here referred to a single simulated pedigree (the empirical one). The empirical results were usually within the 95% of the distribution of simulation results across replicates, concurring with some of their major findings. For example, estimates of the average inbreeding coefficient based on F VR2 and F LH 2 when the allele frequencies in the base population were known ( F , yellow diamonds) had the largest differences with F IBD and the largest variances in F values (Fig. 3a, c). In addition, estimates from these estimators were highly negatively correlated (Fig. 3g), as predicted by the simulations. When base population frequencies were assumed to be equal to 0.5 ( F q05 ), substantially overestimated ΔID based on F IBD (Fig. 3e) were obtained, in agreement with the overestimations found in the scenario EC in the simulations (Fig. 1, Additional file 5: Fig.  S3 and Additional file 6: Fig. S4).

Comparison between simulated and empirical results in an Iberian pig population
Because ROH estimates of inbreeding consider more ancient inbreeding (i.e., inbreeding relative to a more ancient generation) than that of cohort 17, Fig. 3b, d, f, h shows simulation and observed results assuming that the base population is the initial cohort 0. The average inbreeding coefficient based on F ROH −5 was very close to that based on F IBD , while that based on F ROH −1 was higher, as it may also consider more ancient inbreeding. This is well illustrated in Fig. 4, which shows the simulated and observed average inbreeding coefficients based on F IBD , F PED , F ROH −1 , and F ROH −5 across generations. Estimates of inbreeding based on F ROH −5 were very close to F IBD and F PED values. For the simulations, the average number of ROH per individual considering all individuals between cohorts 17 and 23 was about 86, with an average length of 10 Mb for ROH > 1 Mb, and about 60 with an average length of 13 Mb for ROH > 5 Mb. These values are similar to those found in the real data: an average of 77.7 ROH per individual with a mean length of 9.99 Mb for ROH > 1 Mb, and 47.9 ROH with a mean length of 14.04 Mb for ROH > 5 Mb. There was, however, some discrepancy between simulated and empirical estimates of ΔID based on F ROH , with the empirical estimates of the latter being higher (Fig. 4f ).

Discussion
Proper estimation of the inbreeding coefficient of individuals is necessary to estimate inbreeding depression and effective population size in domestic and wild populations. The possibility of obtaining data from large numbers of SNPs for an increasing number of species allows the estimation of genomic inbreeding in the absence of pedigree records or as a complement to them. Here, we carried out computer simulations to compare some of the most popular marker-based estimators of F that are readily obtained from widely-used software such as PLINK and GCTA under different scenarios, including equalization of parental contributions, and artificial and natural selection. Our results show that, all estimators generally had high correlations with F IBD and precise estimates of ΔID when the allele frequencies in the base population are known, with F VR2 and F LH 2 having the lowest performance (Fig. 2). When allele frequencies in the base population are not known, F LH1 generally had the largest correlation with F IBD , and F YA2 generally gave the most accurate estimates of ΔID. The F VR2 and F LH 2 estimators had very low correlations with F IBD and low accuracies to estimate ΔID in most scenarios. The assumption of allele frequencies of 0.5 in the calculation of the estimators ( F q05 ) resulted in high correlations with F IBD and relatively good estimates of ΔID, except in the EC scenario, for which a substantial overestimation was found. Estimates of ΔID obtained with simple homozygosity measures ( F HOM ) resulted in large overestimations of ΔID, although F HOM has a correlation of 1 with F LH 1 ; however, its variance is much lower than that of F IBD (Fig. 1), resulting in estimates of ΔID based on F HOM to be biased upwards. This agrees with previous studies that showed that estimates of ΔID based on homozygosity measures of F are generally larger than those based on other estimators of F [37,49]. Estimates of inbreeding based on ROH gave precise estimates of F IBD and ΔID, as also shown by previous studies [9,[27][28][29][30]47]. Finally, estimates of F based on pedigrees ( F PED ) were less correlated with F IBD than molecular estimates, but provided estimates of ΔID with little bias. Our results for the pedigreed pig population showed generally good agreement between simulated and empirical results, supporting the above findings.

Estimation of the inbreeding coefficient
When allele frequencies for a base population of reference that consists of noninbred and unrelated individuals are known, all marker frequency-based estimators are expected to give almost unbiased estimates of F IBD [19,22,31,32], as observed in our simulation results (Fig. 1). The slightly lower average inbreeding coefficients based on molecular estimates compared to F IBD (Fig. 1) is explained by the negative excess of heterozygotes expected for finite population sizes, which amounts approximately to − 1/2N [69,70]. Note that, as the use of genomics becomes more common, it will be more frequent to have several previous generations genotyped and, therefore, to have information available for a base initial population of reference. Allele frequencies for a base population can also be estimated. Gengler et al. [33] and, more recently, Aldrige et al. [71] have proposed practical methods to calculate gene content using a linear regression. However, these methods are not feasible for deep pedigrees, and it is standard practice to use allele frequencies for the current genotyped population because of the ease of computation. Even if base population frequencies are available, the correlation with F IBD is not the same for all estimators. The F LH 1 , F VR1 , and F YA1−2 estimators provided accurate predictions, but F VR2 and F LH 2 did not. The latter (called F II in the software PLINK and GCTA) is almost never used in practice. However, F VR2 ( F I in PLINK and GCTA) has been considered on many occasions [12,41,43,44,50,58]. Our results argue against the use of both these estimators, as they showed the lowest correlations with F IBD and the most biased estimates of ΔID in the majority of scenarios, as summarized in Fig. 2. When allele frequencies in the base population are not known, and those of the current population are used instead, which is the most frequent situation, the most accurate marker-frequency-based estimator of F IBD was F LH 1 (or F HOM or F q05 , which both had a correlation of 1 with F LH 1 ). F LH 1 had the highest correlation with total homozygosity [32], and it is also the most useful to indicate if genetic diversity has been lost or gained in specific regions of the genome from the past generations to the current one [12]. For example, an increase in the value of local F LH 1 indicates a loss in heterozygosity, while a decrease indicates a gain in heterozygosity. In contrast, F VR1 , F VR2 , and F YA2 do not provide useful information of whether heterozygosity has declined or increased in genomic regions [12].
Estimates of the inbreeding coefficient based on ROH were very accurate estimators of F IBD , although they depended on the length of the fragments considered. For the pig data, Fig. 4 shows that estimates of F ROH −5 were very close to the expected values of F IBD and F PED relative to the initial cohort, while F ROH −1 provide larger estimates of inbreeding. Because the average rate of recombination for pigs is about 0.8 cM per Mb, F ROH −5 values refer to events of inbreeding that occurred since about 12.5 generations ago [9], which probably includes most of the inbreeding that occurred since the foundation of the strain. However, F ROH −1 refers to events of inbreeding that occurred further in the past (on average about 62.5 generations), which includes inbreeding that occurred before foundation of the strain, and explains why the mean F ROH −1 in cohort 0 was greater than 0 (0.1).
Estimates of inbreeding based on current allele frequencies had very variable correlations with each other and with pedigree-based estimates of inbreeding. For example, the non-exhaustive compilation of correlations in Table 2 shows that the correlation of the estimates of inbreeding based on homozygosity ( F LH 1 or F HOM ) with F PED ranged from 0.1 to 0.9, with an average of about 0.5. Correlations of other marker-based estimators with F PED and with each other were also rather variable, with some negative values, usually involving F VR2 . These empirical results from Table 2 agree with our results, which also showed the lowest correlation between F PED and F VR2 , both in the simulations [see Additional file 9: Fig. S6] and in the pig data (Fig. 3). Correlations between molecular estimators of F were generally high, except for some correlations that involved F VR2 particularly with F LH 1 and F LH 2 , which were often low or negative ( Table 2). This is also evident from our simulations results and the pig data ([see Additional file 9: Fig. S6] and Fig. 3), as well as from simulation results by Zhang et al. [34], where correlations of F VR2 with F PED , F LH 1 , and F YA1 were negative. Moreover, minima and maxima of individual inbreeding coefficients were within the expected ranges in our study [see Additional file 7: Table S3]), except for some values higher than 1 or lower than -1 for F VR2 and F LH 2 , which are the estimators with the poorest performances. Therefore, the different estimators allow for a wide range of possible values of F, which likely depends on many factors such as the species, mating system, family structure, density and quality of markers, population size, etc. In addition, empirical results suggest that genotype imputation largely increases the variance of the estimators, thus producing substantial overestimates of average inbreeding coefficients and individual values well outside the expected ranges [44].

Estimation of the rate of inbreeding depression
Regarding estimation of ΔID when current allele frequencies are used, the most accurate molecular marker frequency-based estimator of inbreeding appears to be F YA2 . This is the estimator that results in a variance of F values that is closest to that of F IBD (see VF in Fig. 1, Additional file 5: Fig. S3 and Additional file 6: Fig. S4), and the lowest variance of individual estimates among all molecular estimators, as deduced theoretically by Yang et al. [36] and confirmed by other simulation [32] and empirical [42] results. F YA2 was also the estimator with the greatest power to detect inbreeding depression [see Additional file 8: Fig. S5]. There has been a number of simulation studies that have investigated the performance of different measures of inbreeding to estimate ΔID when allele frequencies of the current population are assumed [32,46,47]. These studies suggest that F YA2 and F ROH generally provide good estimations of ΔID, depending on the population size, while F LH 1 slightly underestimated ΔID and F VR2 severely underestimated ΔID. In fact, F YA2 and F ROH were highly correlated ( Table 2) and previous work has shown that F YA2 and F ROH have approximately the same power to detect inbreeding depression [13,37], as also confirmed by our simulation results [see Additional file 8: Fig. S5]. It has also been found [32] that F YA2 has the highest correlation with the phenotypic values of individuals and with the homozygous mutation load [9], as a proxy for fitness. However, results from previous simulations demonstrate that F YA2 can substantially overestimate ΔID for relatively small population sizes [32,47]. For example, an overestimation by about 28% was observed for a population of size of N = 500 [32]. This substantial overestimation of ΔID based on F YA2 was not observed in our results, where F YA2 was always rather accurate. The reason for this discrepancy is that the previous studies analysed samples that were taken directly from a relatively large population, while we simulated a number of generations with a much smaller N (20 or 100) prior to the generation that was analysed, and mutation was not considered in the simulations. Thus, in our simulations, many rare alleles were likely lost by genetic drift during the period of reduced census size, which may have improved the estimation of the rate of inbreeding Table 2 Compilation from the Literature of correlations between individual measures of inbreeding from empirical data based on marker-based estimators of inbreeding assuming allele frequencies in the current generation and based on runs of homozygosity (ROH) The compilation is not intended to be exhaustive, but to give a general picture of the empirical correlations found in different studies with several marker-based estimators of inbreeding. F VR1 and F VR2 : Estimators described by VanRaden [19]. F YA1 and F YA2 : Estimator proposed by Yang et al. [20] and modified by Zhang et al. [34]. F LH1 and F LH2 : Estimators derived by Li and Horvitz [16]. − 0.09 [39,41,52,54,55] 0.42 [13,39,42,49,55,57] 0.50 [13,14,39,41,42,49,51,52] 0.32 [55] 0.50 [13,14,39,41,49,51] 0. 10 [55] depression. In order to check this, we also analysed the case with N = 100 with the RC scheme at generation 0, i.e., immediately after sampling from the large original populations with N = 500. Resulting estimates of ΔID are given in Additional file 11: Fig. S8, which show that, with no minor allele frequency (MAF) filter (red bars), F YA2 overestimated ΔID by about 28%, F LH 1 underestimated ΔID by about 30%, while F ROH produced almost unbiased estimates and F VR2 and F LH 2 gave large underestimations. These results are in full agreement with those previously reported [32]. However, if these analyses are done when rare alleles were removed by applying a MAF filter of 0.01 or 0.05 to the data, overestimation of ΔID by F YA2 was reduced to 3% for the 0.01 filter and became an underestimation by 15% for the 0.05 filter [see Additional file 11: Fig. S8]. Thus, part of the overestimation of ΔID that occurred with F YA2 can be due to the contribution of SNPs with very low frequencies. Estimates of ΔID based on F VR1 , F YA1 , and F LH 1 , i.e. those based on ratios of averages over loci, were not affected by applying an MAF filter, while those based on F VR2 , F LH 2 , and F ROH were increased by an increase of the MAF threshold [see Additional file 11: Fig. S8].
The results of our study are concordant with those of Wang [23]. Under the conditions considered in our simulations (genome size of 18 Morgans, number of SNPs > 30,000, population sizes < 100, and number of generations > 20), molecular estimates were more precise estimates of F IBD and ΔID than those based on pedigree. Wang [23] considered random mating populations under a model of deleterious mutations on fitness, where the molecular inbreeding metric ( F LH 1 ) was calculated using current allele frequencies. We reached similar conclusions (Fig. 2a) for other molecular estimators of inbreeding and scenarios with equalization of contributions and artificial selection. Our results also concur with those of Forutan et al. [48], who compared estimates of inbreeding obtained from pedigree and from molecular data, i.e., F VR1 and F ROH , considering random mating and artificially selected populations, and a base population simulated assuming allele frequencies equal to 0.5 or uniformly distributed between 0 and 0.5. They showed that pedigree-based estimates of inbreeding underestimate inbreeding in selected populations, a result that we also observed in our simulations [see Additional file 6: Fig.  S4]. This is as expected because selected individuals tend to be more related than the mean of the population [72]. They also found that F VR1 estimates were close to F IBD estimates if the allele frequencies of the base population are known, but assuming a constant allele frequency of 0.5 resulted in overestimation when using F VR1 , although the correlation of these estimates with F IBD was close to 1. In fact, high correlations (> 0.6) between inbreeding values obtained based on F PED and F q05 have been found in several empirical studies [35,37,48]. This agrees with our results, which also showed that F q05 generally provided accurate estimates of ΔID (Fig. 2b), except for the EC scheme, where it resulted in overestimation (Fig. 1, Additional file 5: Fig. S3 and Additional file 6: Fig. S4), which was also reflected by the pig data (Fig. 3), a population that was intended to follow an EC protocol.
Pedigree-based inbreeding ( F PED ) showed rather accurate mean estimates of F IBD (see Fig. 1, Additional file 5: Fig. S3 and Additional file 6: Fig. S4), except in the case of artificial selection (see Additional file 6: Fig. S4 and [48]), but had only a moderate correlation with F IBD (Fig. 2a), as also found by Wang [23]. Estimates of ΔID by F PED were on the whole rather accurate (Fig. 2b), which concurs with the meta-analysis of Doekes et al. [73], who found moderate to high correlations of the estimates of ΔID based on F PED with those obtained based on molecular estimators of inbreeding: i.e., about 0.5-0.6 for F ROH , about 0.6-0.7 for F VR2 or F YA2 (increasing to 0.8-0.9 for F YA2 when base population frequencies of 0.5 were used), and about 0.6-0.7 for F LH 1 . According to Wang [23], molecular estimates of inbreeding are more accurate predictors of F IBD than pedigree-based estimates when population sizes are smaller than ~ 200 (see his Fig. 3). This agrees with other simulation results assuming small population sizes [28]. However, Wang [23] also showed that for larger population sizes (> 200), the accuracy of pedigree-based estimates of inbreeding increased and became higher than that from molecular estimates of inbreeding. Whereas the accuracy of the molecular estimates of inbreeding was reduced with increasing population size, that from pedigree-based estimates remained invariable with increasing population size (see his Fig. 3). In addition, power to detect inbreeding depression became greater for pedigree-based estimates than for molecular estimates when population sizes were larger than about 500 (see Figure 5 of Wang [23]). Therefore, for large population sizes, under some situations, estimates of inbreeding and ΔID obtained from pedigrees can be as reliable, or even more reliable, than estimates based on molecular markers, although pedigrees for such large populations may only be available in domestic species. In most situations, therefore, estimates of inbreeding and ΔID should be more reliably obtained from molecular markers than from pedigrees.
Our simulations implied that the effective population size ( N e ) was close to the census size of breeders ( N ), except for the EC scenario, for which N e ≈ 2N ; see e.g., [4], p. 110). Thus, the results referred to scenarios with relatively small N e (≤ 100), with average inbreeding coefficients reaching relatively high values ( F ≈ 0.1 − 0.4). This is a common situation in many animal breeding populations [74] and populations of conservation concern [75]. Other simulation studies have focused on populations with much larger population sizes (e.g., N > 10,000 individuals) and these have shown that molecular estimators of inbreeding perform also rather well [23,32,46].
In order to estimate inbreeding depression for fitness, we assumed a model of partially recessive deleterious mutations with average mutational parameters based on empirical data (see, e.g., [4]; p. 152-161). We also considered an alternative model with a larger number of mutations with much lower effects (Table 1 and Additional file 1: Fig. S1) but this did not affect the main findings of the study. However, we did not consider models of overdominance or epistasis. Overdominance can make some contribution to inbreeding depression but this is generally considered to be minor compared to partial dominance [76][77][78]. Nevertheless, further studies should be devoted to the investigation of the performance of the molecular estimators of inbreeding under other genetic models.

Conclusions
If the allele frequencies of the base population are known, all marker frequency-based estimators of inbreeding and inbreeding depression are reasonably highly correlated with F IBD and provide accurate estimates of ΔID, except for F VR2 and F LH 2 . If the allele frequencies of the base population are not known, F LH 1 is generally the estimator that is best correlated with F IBD and F YA2 is the best estimator of ΔID. F ROH is a very accurate estimator of F IBD in most situations, while F VR2 and F LH 2 perform very poorly in almost all scenarios. Estimates that are obtained assuming a constant frequency of 0.5 ( F q05 ) give highly biased estimates of F IBD but are highly correlated with F IBD and provide reasonably good estimates of ΔID. Estimates from simple frequencies of homozygous markers ( F HOM ) cannot be used to estimate ΔID. p k q k Because, p 2 kg1 = p kg1 and p 2 kg2 = p kg2 , then, Equation (3), which was also shown in the Appendix of Caballero et al. [32], indicates that the estimator is based on the squared deviation of the allele frequencies of individuals from the population mean minus the expected variance within individuals. This contrasts with the estimator of VanRaden [19], which is based on the correlation of allele frequencies between individuals and does not include the correction for the expected variance within individuals in the numerator of its expression, (see Appendix of Caballero et al. [32]). After some algebra on Eq. (3), (2p ki − p k ) (2p ki − p k ) − p kg1 + p kg2 + p kg1 p kg2 p k q k = 1 S S k=1 −(2p ki − p k )p k + p kg1 p kg2 p k q k .