An empirical evaluation of the estimation of inbreeding depression from molecular markers under suboptimal conditions

Abstract Inbreeding depression (ID), the reduction in fitness due to inbreeding, is typically measured by the regression of the phenotypic values of individuals for a particular trait on their corresponding inbreeding coefficients (F). While genealogical records can provide these coefficients, they may be unavailable or incomplete, making molecular markers a useful alternative. The power to detect ID and its accuracy depend on the variation of F values of individuals, the sample sizes available, and the accuracy in the estimation of individual fitness traits and F values. In this study, we used Drosophila melanogaster to evaluate the effectiveness of molecular markers in estimating ID under suboptimal conditions. We generated two sets of 100 pairs of unrelated individuals from a large panmictic population and mated them for two generations to produce non‐inbred and unrelated individuals (F = 0) and inbred individuals (full‐sib progeny; F = 0.25). Using these expected genealogical F values, we calculated inbreeding depression for two fitness‐related traits, pupae productivity and competitive fitness. We then sequenced the males from 17 non‐inbred pairs and 17 inbred pairs to obtain their genomic inbreeding coefficients and estimate ID for the two traits. The scenario assumed was rather restrictive in terms of estimation of ID because: (1) the individuals belonged to the same generation of a large panmictic population, leading to low variation in individual F coefficients; (2) the sample sizes were small; and (3) the traits measured depended on both males and females while only males were sequenced. Despite the challenging conditions of our study, we found that molecular markers provided estimates of ID that were comparable to those obtained from simple pedigree estimations with larger sample sizes. The results therefore suggest that genomic measures of inbreeding are useful to provide estimates of inbreeding depression even under very challenging scenarios.


| INTRODUC TI ON
Wild and domestic small populations are often affected by inbreeding (see, e.g., Hasselgren & Norén, 2019), which can lead to a decline in the average fitness-related traits, known as inbreeding depression (Doekes et al., 2021;Hedrick & Garcia-Dorado, 2016;Lynch & Walsh, 1998, Chap. 10). Inbreeding depression (ID) is of paramount relevance for the resilience of populations and the success of breeding programs (Frankham, 2003(Frankham, , 2005Robertson, 1961;Wright et al., 2008), and occurs primarily due to the increased homozygosity of (partially) recessive deleterious mutations that are concealed in heterozygosis in a non-inbred population. Inbreeding depression can be estimated by calculating the slope of the linear regression of the phenotypic values of individuals on their inbreeding coefficients (Lynch & Walsh, 1998). These coefficients can be obtained from pedigree data (Wright, 1969). However, pedigrees are not available for many organisms, particularly in wild populations (Pemberton, 2008), and molecular markers can be used instead. With the increasing availability of a large number of molecular markers, such as single nucleotide polymorphisms (SNPs), it is now possible to estimate the genomic inbreeding coefficient of individuals without requiring knowledge of their pedigree-based relationships (Goudet et al., 2018;Howard et al., 2017;Wang, 2014;Yengo et al., 2017).
In addition, genomic inbreeding has the advantage of providing a more precise measure of the relationships between relatives compared to the pedigree coefficient, as the latter gives just expected values (Howard et al., 2017;Kardos et al., 2016;Wang, 2016).
There are several measures of the inbreeding coefficient of individuals from genomic data. These measures can be obtained either on a SNP-by-SNP basis, which depend on allele frequencies (e.g., Li & Horvitz, 1953;VanRaden, 2008;Yang et al., 2010), or from runs of homozygosity (ROH), that is, regions of the genome that are homozygous for all or the vast majority of nucleotide bases (Broman & Weber, 1999;McQuillan et al., 2008). Numerous simulation and empirical studies have compared the different measures of inbreeding based on markers, showing that some of them are very appropriate estimators of inbreeding depression (Caballero et al., , 2022Keller et al., 2011;Nietlisbach et al., 2019;Wang, 2014;Yengo et al., 2017).
Both empirical and simulation studies indicate that estimates using data from different and relatively distant generations, implying a substantial variation in F values, increase the power to estimate ID (e.g., Bérénos et al., 2016;Doekes et al., 2019;Martikainen et al., 2017;Saura et al., 2015;Silió et al., 2013;Sumreddee et al., 2019). However, in many cases, historical data may not be available, and the F estimate is obtained using individuals from the contemporaneous population. In such situations, the variation in F values among individuals can be low, particularly for large sized populations, which may compromise the estimation of ID. Furthermore, the power to detect ID strongly depends on the sample sizes (Wang, 2014) and the reliability of the estimates of individual fitness, which can be subject to substantial error and sources of environmental variation due to their low heritabilities (see, e.g., Caballero, 2020, p. 57).
In this study, we aimed to evaluate inbreeding depression for two fitness traits of Drosophila melanogaster under restricted experimental conditions using individual F values estimated from molecular markers. Firstly, individuals were sampled from a large panmictic population at a single generation, which would usually lead to low variation in F values. Secondly, the sample size for analysis was rather limited with only 17 pairs of individuals. Lastly, the traits measured depended on both males and females, but only males were sequenced and had their inbreeding coefficients estimated. Our results show the usefulness of genomic measures of inbreeding for estimating inbreeding depression, despite the strong experimental limitations.

| ME THODOLOGY
In brief, we estimated the inbreeding depression rate for two fitness traits (pupae productivity and competitive fitness) in a large panmictic laboratory population of D. melanogaster, following a fullsib mating design. One hundred non-inbred pairs of individuals and 100 inbred pairs (full-sib progeny) were established and scored for the traits when possible. This allowed an estimate of the inbreeding depression to be obtained from the expected pedigree estimates of inbreeding (0 and 0.25 respectively). DNA was individually extracted for a set of 17 males from the non-inbred pairs, and 17 males from the inbred pairs, and whole-genome sequencing was carried out obtaining a large number of SNPs. Estimates of the inbreeding depression rate were then obtained from the linear regression of the individual fitness values on the genomic measures of inbreeding obtained by means of two estimators. The estimated ID rates obtained from the markers were compared with those obtained from the simple pedigree values of expected inbreeding considering a much larger sample size.

| Base population, breeding design and fitness traits evaluated
A population of D. melanogaster was founded in October 2018, derived from a wild population located in a wine cellar in Ponteareas (Pontevedra, Spain). This population was maintained in the laboratory for 50 generations at a constant population size of around N ≈ 2500 individuals, distributed in 32 bottles with ~80 individuals per bottle. The laboratory conditions were kept invariant over generations, with a constant temperature of 25°C and continuous lighting. To ensure a near panmictic scenario, circular mixing of individuals from consecutive bottles was carried out each generation.
Individuals introduced in each bottle were removed after a week in order to avoid the overlap of generations. A summary of the overall experimental design is presented in Figure 1.
An evaluation of fitness and inbreeding depression was carried out at generation 50, by sampling at least three virgin females and three males from each bottle and randomly placing pairs of individuals in single vials, obtaining a total of 100 breeding pairs with the subsequent separate offspring. This was considered a transition generation to exclude maternal effects (generation t = 0). Two separate schemes were derived from the offspring: First, a non-inbred scheme, where individuals were mated randomly in individual vials for two generations (100 breeding pairs) but avoiding inbred matings. This was achieved by mating individuals from distant vials. For example, at generation t = 1, a female from vial 1 was mated with a male from vial 48, a female from vial 2 with a male from vial 49, and so on, following a circular scheme (see Figure 1). Analogously, at generation t = 2, a female from vial 1 was mated with a male from vial 20, a female from vial 2 with a male from vial 21, and so on. Thus, adult individuals at generation t = 2 were expected to be unrelated with an expected inbreeding coefficient F PED = 0, relative to the base population, in this case the generation 50 of the large population from which the experiment started. Second, an inbred scheme was derived by mating full-sib individuals for two generations. Thus, adult individuals had an expected inbreeding coefficient F PED = 0.25 at generation t = 2, relative to the base population, but individuals from different families were also expected to be unrelated among each other.
Two fitness-related traits were simultaneously evaluated at generation t = 2 for the two sets of individuals analysed. On the one hand, pupae productivity (P) was obtained as the mean number of pupae produced from each breeding pair in each scheme after 11 days from mating (t = 2). On the other hand, to evaluate competitive fitness (hereafter referred simply to as fitness, W), females were taken from the vials of the two schemes 3 days after the mating (t = 2) and were placed together with Curly-type females (previously randomly mated in bottles) in groups of eight females per vial, that is, four independent wild-type females (i.e., from four different vials) plus four Curly-type females per vial (see Figure 1). In the inbred scheme (F PED = 0.25), it was not possible to include fully unrelated wild-type females in all evaluation vials, and about half of the vials included two wild-type females which were cousins, and therefore had an expected coancestry coefficient of f = 0.25. After 14 days, wild-type offspring and Curly-type offspring were counted from each evaluation vial and fitness (W) was measured as the ratio Wild/(Curly + 1). Note that both fitness measures are complex traits involving mating success, parental fecundity and egg-to-pupae viability in the case of productivity (P), and also egg-to-adult viability including competition against the Curly strain in the case of fitness (W; see, e.g., López-Cortegano et al., 2016). Thus, both traits measured in the progeny depend on both male and female parents. In the case of pupae productivity, it was previously shown that about half the ID shown in the trait is due to the female parents (Avila et al., 2013). F I G U R E 1 Experimental design followed to evaluate fitness and inbreeding depression. Two schemes were derived from the base population to achieve an expected inbreeding coefficient of F PED = 0 and F PED = 0.25. Fitness was measured as mean pupae productivity (P for both the non-inbred and inbred schemes) produced after 11 days since the last mating (t = 2), and as competitive fitness (W for both the non-inbred and inbred schemes, as the ratio of Wild-type to Curly-type offspring) 14 days after placing those mated females at t = 2 in the presence of curly females in groups of four (i.e., four Wild-type females plus four Curly-type females). Sequenced males are shown in red (see main text for details on how males were selected for sequencing).
The male parents from about 60 randomly chosen non-inbred pairs and about 60 inbred pairs at generation t = 2 were taken to be frozen for later DNA extraction and sequencing (see next section), but only a subset of 17 non-inbred and 17 inbred males were sequenced.
Female parents were not considered for sequencing to avoid offspring contamination. Note that the individual pair productivity (P) corresponds to the pair, that is, to the male (sequenced) and the fe-

| DNA extraction, sequencing and SNP calling
As mentioned above, breeding males from generation t = 2 were removed from their corresponding vials 3 days after mating, frozen with liquid nitrogen and individually stored at −80°C until DNA extraction. DNA extraction was carried out using the Gentra Puregene Cell Kit (Qiagen) with some modifications, including an RNase step.
From all males frozen, we randomly chose 17 from the non-inbred scheme (F PED = 0) and 17 from the inbred scheme (F PED = 0.25) from those meeting the following requirements: (i) a minimum DNA concentration of 10 ng/μL (in a final volume of 10 μL), measured in a Qubit fluorimeter (Thermo Fisher Scientific), (ii) a minimum DNA integrity number (DIN) of 8 measured in an Agilent TapeStation system, and (iii) the presence of a phenotypic value for both pupae productivity (P; those individuals showing a phenotypic value equal to zero were excluded) and fitness (W; avoiding males from the same group, that is, only one male was selected per group of four families as indicated in Figure 1, where the males that were sequenced are represented in red). Nextera XT DNA library preparation and genome sequencing were carried out at the NimGenetics Genomics Service, preparing a single Illumina library for each individual fly. A total of 34 DNA samples (17 from individuals with F PED = 0 and 17 from individuals with F PED = 0.25) were used for 2 × 150 bp paired-end sequencing on the Illumina Novaseq 6000 instrument to a mean read depth of ~50 times.
Paired-end reads were processed from the FASTQ files to finally obtain the Variant Call Format (VCF) files with data for filtered biallelic autosomic SNPs. Briefly, adapters were removed with Trimmomatic (Bolger et al., 2014) using the Nextera adapter list, and reads were filtered according to their quality and size with the ERNE-FILTER tool (default options except for --min-size = 36; Del Fabbro et al., 2013). The resulting reads were mapped to the v. 6.14 D. melanogaster reference genome from Flybase (www.flyba se.org/) using the BWA-MEM algorithm from the Burrows-Wheeler Alignment (BWA) software (Li, 2013). After mapping, PCR duplicates were removed and the alignments filtered for a minimum map quality score of 20 with SAMtools (Danecek et al., 2021). Potential variant sites were called with the HaplotypeCaller tool from the Genome Analysis Toolkit (GATK; Van der Auwera & O'Connor, 2020). The resulting GVCF (Genomic Variant Call Format) files were combined into one using the GATK CombineGVCFs tool. The GATK GenotypeGVCFs tool was used to obtain genotype information for all variants.
Multiallelic and monomorphic SNPs, as well as indels and those SNPs within 10 bp of an indel, were removed with BCFtools (Danecek et al., 2021). Highly repeated sites and low information areas were also excluded with VCFtools (Danecek et al., 2011). Around 1,140,000 SNPs from autosomal chromosomes were finally available for inbreeding analysis. Detailed information about the sequencing quality parameters, such as the mean mapping quality (58.5 on average), mean coverage (49× on average), percentage of the genome covered by at least one read (93.5% on average) and total number of mapped bases (5.4 × 10 9 on average) for each individual can be found in Supplemental File S1, along with DNA extraction data.

| Estimation of inbreeding depression from molecular inbreeding measures
Two estimators of the inbreeding coefficient were obtained with PLINK (version 1.9) from the SNPs available for individuals with F PED = 0 and for individuals with F PED = 0.25 separately. The first estimator is a SNP-by-SNP estimator obtained with the command -ibc (equivalent to the -ibc from GCTA; Yang et al., 2011) and calculated as where S is the total number of markers, x k the number of minor alleles of marker k (0, 1 or 2 copies) and p k the current frequency of the minor allele of marker k (Yang et al., 2010). F YAN (called F III in PLINK) is based on the correlation between uniting gametes.
Estimates of F were also obtained from runs of homozygosity (ROH), calculated as where ΣL ROH is the sum of the lengths of all ROHs (above a specified minimum length) that cover the genome of an individual, and L auto the length of the autosomal genome covered by SNPs (McQuillan et al., 2008), that is, 119 Mb for autosomal chromosomes considering that around 90% of the reference genome was sequenced. First, LDbased pruning was carried out to remove highly linked SNPs (r 2 > 0.9) using the --indep-pairwise 50 5 0.9 PLINK option, as recommended by Howrigan et al. (2011). Around 420,000 SNPs remained after pruning. Then, F ROH was obtained with the --homozyg command using the default options to define a ROH, among which stands out a minimum number of SNPs of 100, a minimum density of 1 SNP per 50 Kb and a scanning window of 50 SNPs (with 1 heterozygote allowed per window). Two minimum ROH lengths were applied, 0.1 Mb (F ROH-0.1 ) and 1 Mb (F ROH-1 ).
Among the various available measures for estimating inbreeding, these two molecular measures have been shown to be particularly reliable for estimating the rate of ID (

| Expectation of the rate of inbreeding depression estimated by the molecular inbreeding measures given the experimental design
Because F estimates from molecular markers were only obtained for males, while the phenotypic values correspond to the pair, that is, the male and the female, in the productivity measure, and to four pairs in the case of fitness, we might expect a partial estimate of the pedigree ID. Following the results from Avila et al. (2013), we assume that both traits depend equally on both parents. In that case, estimates from unrelated individuals (F PED = 0) are expected to be one half the value estimated if both sexes were considered in the case of productivity, and one eight in the case of fitness. Analogously, estimates from full-sib individuals (F PED = 0.25) are expected to be a fraction around 7/8 and 9/32 for productivity and fitness, respectively (see Supplemental File S2).

| RE SULTS
Mean molecular measures of the inbreeding coefficient are shown in Table 1. SNP-by-SNP estimates (F YAN ) provided values close to zero and slightly negative for individuals with F PED = 0, as expected.
ROH-based estimates provided larger values, especially when a minimum length of 0.1 Mb (F ROH-0.1 ) was applied, compared to a minimum length of 1 Mb (F ROH-1 ) and, to a larger extent for individuals with F PED = 0.25, as expected. Standard errors (SEs) were also larger for shorter ROH lengths. For all estimates, standard deviations between F values of individuals were larger for F PED = 0.25 than for F PED = 0. Table 2 gives the information of the ROH fragments found.
Individuals with F PED = 0.25 presented a larger proportion of the genome with ROH, with 5% more fragments per individual on average than individuals with F PED = 0 for a minimum ROH length of 0.1 Mb, and about twice the number of fragments when a minimum ROH length of 1 Mb was applied. No differences in SNP density within ROH were observed between ROH lengths and origin of individuals.
Pearson correlations between F estimates (Table 3; p-values in Table S1, Supplemental File S3) were very high in all cases (>0.70), but to a greater extent for individuals with F PED = 0.25. Table 4 presents the Pearson correlations between molecular F values and the two fitness traits, Log(P) and Log(W). Overall, F YAN showed the highest (negative) correlation with Log(P) and Log(W). Molecular F estimates were more correlated with Log(P) than with Log(W), at least for individuals with F PED = 0, as expected given the experimental design, although the opposite occurred with individuals with F PED = 0.25. Despite these trends, correlations were not significant in any case (p-value >0.08; Table S2 in Supplemental File S3).
Inbreeding depression (ID) was obtained as the slope of the linear regression of the phenotypic values, Log(P) and Log(W), on the inbreeding coefficient estimates. All estimates of the rate of inbreeding depression are summarized in Table 5. Regarding Log(P)  (Table S3). Although there were some differences between methods, there was a general concordance between them. When the estimates of ID were obtained grouping the two sets of data (i.e., with 34 samples instead of 17), these were in general similar to those obtained for individuals with F PED = 0.25 (black regression lines in Figures 2 and 3), with the difference that now all ID estimates were significantly different from 0 irrespective of the test method (p-value <0.041; Table S3 in Supplemental File S3), and were closer to the ID estimates obtained from F PED measures in the case of productivity ( Table 5).
The above results were obtained considering whole genome sequencing data. In order to assess the performance of the estimators when a lower number of SNPs is available, mean F values and ID estimates were obtained from a random subset of SNPs (from 500,000 down to 5000; Tables S4 and S5 in Supplemental File S3).
Mean F values were not particularly affected by the number of SNPs, except for ROH-based estimates (

TA B L E 4
Pearson correlation between different molecular measures of the inbreeding coefficient (F; see text) and the phenotypic values for two fitness traits: pupae productivity (P) and fitness (W), for the two sets of individuals (non-inbred, F PED = 0; and inbred, F PED = 0.25). while F ROH-0.1 were slightly decreased. For a number of SNPs below 20,000, both estimates offered virtually the same values, as short ROH were more difficult to detect (e.g., for F ROH-0.1 with 5000 SNPs, the minimum length of a segment detected was around 2 Mb; not shown). Analogously, ID estimates were not particularly affected by the number of SNPs when using SNP-by-SNP estimators (Table S5) Table 5.
−0.355 to −0.034. We also did not find large differences in ROHbased IDs when modifying the default parameters in PLINK (Table   S6 in Supplemental File S3).

| DISCUSS ION
The power of molecular markers to estimate the rate of inbreeding depression is influenced by several factors, such as the level of variation in individual F values, the accuracy of F estimation, the sample sizes available, and the reliability of fitness measures. Empirical studies on the performance of molecular markers in estimating ID typically involve samples from multiple generations, which allows for a wide range of inbreeding coefficients among individuals, which facilitates estimation accuracy. For example, most studies include pedigrees spanning seven to 12 pedigree generations (Antonios et al., 2021;Bérénos et al., 2016;Doekes et al., 2019;Ferenčaković et al., 2017;Makanjuola et al., 2020;Martikainen et al., 2017;Zhang et al., 2022), with some studies covering more than 20 generations (Saura et al., 2015;Silió et al., 2013;Sumreddee et al., 2019) and a few spanning fewer than five generations (Hidalgo et al., 2021;Pryce et al., 2014). These studies are also characterized by involving large sample sizes, exceeding 1000 individuals and sometimes reaching over 20,000 individuals (Doekes et al., 2019;Makanjuola et al., 2020). While results vary due to species differences, trait variation, and study-specific factors, molecular estimators generally  Table 5.
such extreme scenarios. Although many of the estimates were nonsignificant with the small sample sizes considered, they all showed a trend compatible with inbreeding depression and estimated values of ID close to those expected, at least for productivity. All estimates became significant when we pooled the two sets of individuals (inbred and non-inbred) and obtained a higher variation in individual F values. Thus, our study highlights the potential of molecular markers for estimating inbreeding depression in situations with low sample sizes of unrelated individuals obtained from a single generation.
We focused on using molecular markers to estimate inbreeding depression in D. melanogaster, a commonly used model species. We found rates of inbreeding depression of 3.27 for productivity and 3.75 for competitive fitness. The estimate for productivity was obtained from the regression of log productivities on F values in order to compare pedigree and molecular estimates of ID. An alternative calculation, carried out in other studies, is to obtain the estimate from the ratio of the log mean productivities. In this case, ID = log(P-   (Agrawal & Whitlock, 2011;Caballero, 2006Caballero, , 2020Manna et al., 2011), with lower values for strongly deleterious mutations than for milder ones (Agrawal & Whitlock, 2011;Caballero, 2006;Huber et al., 2018;Mukai et al., 1972;Simmons & Crow, 1977). It should be noted that the above estimates refer only to a fraction of mutations of relatively large effect, which can be detected in MA experiments (García-Dorado et al., 2004;Halligan & Keightley, 2009).
This class of mutations, along with that of lethals, which arise at a rate of about 0.015 (Simmons & Crow, 1977), is expected to be the most relevant when it comes to assessing the impact of inbreeding depression and the evolution of fitness in small populations and short periods of time, as shown by comparisons between computer simulations and empirical results (e.g., Caballero et al., 2002;Caballero & Keightley, 1998;Pérez-Pereira et al., 2021). In Drosophila, about half the inbreeding load observed is due to lethal mutations (Lynch & Walsh, 1998, p. 281;Simmons & Crow, 1977). to lack of recombination in males) is a peculiarity that affects the estimation of the rate of ID. It is expected that estimates of inbreeding coefficients will vary due to Mendelian sampling and linkage, with this variation being particularly significant for short genomes (Hill & Weir, 2011). For example, in the case of large genomes such as the human one, the expected standard deviation of the proportion of the autosomal genome shared between full siblings is approximately 0.04 (Hill & Weir, 2011;Visscher et al., 2006). However, this variation is much greater for short genomes (Hill & Weir, 2011 We used two different measures of inbreeding to estimate the inbreeding of individuals in our study. The first measure was the SNPby-SNP estimator of inbreeding (F YAN ), which has shown excellent performance in simulation studies and has been found to be highly accurate in estimating ID in various scenarios (Caballero et al., , 2022Nietlisbach et al., 2019;Yengo et al., 2017). In fact, F YAN had the strongest correlations with both the phenotypes Yengo et al., 2018) and the homozygous mutation load (Alemu et al., 2021;Caballero et al., 2021), which has been suggested to be a proxy for fitness . We also considered estimates obtained from ROH fragments, which have been repeatedly shown to provide good estimates of ID (Caballero et al., , 2022Curik et al., 2014;Keller et al., 2011;Nietlisbach et al., 2019). In comparison, other measures of inbreeding from molecular data, such as those from Li andHorvitz (1953), VanRaden (2008), or the direct SNP homozygosity, are expected to produce more biased estimates of ID than F YAN or F ROH (Caballero et al., , 2022. Average values of the F estimates obtained with molecular data were generally consistent with expectations ( Table 1). SNP-by-SNP-based estimates provided deviations from Hardy-Weinberg proportions or correlations between alleles when the current allele frequencies were used instead of those from the base population, which are usually unknown (Wang, 2014). In random mated finite populations, the expected values are typically slightly negative due to the excess of heterozygotes, approximating −1/(2N) (Kimura & Crow, 1963). Larger positive values are expected due to the loss of heterozygosity as the proportion of full-sib mating increases (Ghai, 1969). ROH-based measures provided larger estimates of inbreeding, as expected, since they are supposed to capture more ancient inbreeding Pryce et al., 2014). We observed a reduction in F estimates (and SEs) for an increase in the minimum length of ROH fragments, with F ROH-1 values close to zero (0.0518 ± 0.0119) for individuals with F PED = 0. Correlations between molecular F estimates were high for all estimators (Table 3), being generally above 0.8, in agreement with those obtained for several domestic populations (see a review of correlations in table 1 of Caballero et al., 2022). Despite the fact that only males were sequenced, and the phenotypic traits depended on both members of the pair (P) or even different pairs (W), we found that F estimates were highly negatively correlated with the phenotypes (pupae productivity and fitness), with correlations around −0.3 ( Table 4).
Estimates of ID were obtained by regressing the phenotypic values on the molecular inbreeding coefficients. When considering a simple design of full-sib mating progeny compared with non-inbred progeny, pedigree-based estimates of inbreeding depression are expected to be accurate. Therefore, we used ID PED (obtained from all individuals with phenotypic values, either sequenced or not) as the standard ID measure against which we compared the molecular estimates. Because the two traits analyzed are complex ones (involving mating success, fecundity, and egg-to-pupae viability for P, plus egg-to-adult viability for W), they are expected to depend on both members of the pairs. As molecular measures of inbreeding were obtained only from males, the estimates of ID from these measures are expected to be 1/2 and 1/8 of those of reference, for P and W, respectively, in the case of non-inbred samples (F PED = 0; Supplemental File S2). For the inbred individuals, because these are the progeny of full sibs, the expected coancestry between the members of the evaluated pair would be f = 0.375, which implies that they share, on average, a 75% of their genome. Thus, the expected fractions of the ID to be estimated by markers was about 7/8 and 9/32 for P and W respectively (Supplemental File S2; Table 5). For productivity, the estimates of ID from F YAN were rather close to the expectations for the set of non-inbred individuals (F PED = 0), whereas both F YAN and F ROH provided good estimates for the set of inbred individuals (F PED = 0.25; Table 5). For fitness, however, the molecular estimates were generally larger than expectations probably because of the large error attached to the estimates.
Estimates of ID were higher for large ROH than for smaller ones (F ROH-1 vs. F ROH-0.1 ), which is consistent with empirical findings. For example, Pryce et al. (2014) reported ID for milk yield for Holstein and Jersey dairy cattle when considering long ROH, but no ID was detected using short ROH. Short ROH have originated from ancient common ancestors (i.e., ancient inbreeding), whose length has been reduced over generations due to recombination, which can lead to overestimations of the current genomic inbreeding. In addition, the effect of genetic purging, which removes deleterious mutations or reduces their frequency, could have mitigated the impact of ancient inbreeding on fitness. As a result, using short ROH may underestimate the extent of current ID. On the other hand, it is important to note that large ROH can also lead to overestimation of ID, as they only indicate recent inbreeding and may therefore underestimate the current genomic inbreeding.
The estimates of ID for fitness (W) were found to be much more variable than those for productivity. This was expected, as the design of the study involved measuring inbreeding in a single male, while four couples were involved in the expression of the trait.
Estimates from inbred individuals, in fact, were substantially larger than expected, and this could be explained by the above reason.
For example, if the family from which the male was sequenced had lower inbreeding and higher fitness than the other three families, the reduced overall fitness observed would be associated with a small F value, which could lead to an overestimation of ID. The sequencing of all individuals in that group could have solved this uncertainty.
An interesting issue arising from Figures 2 and 3, is that, for a given F value, the fitness of inbred individuals (F PED = 0.25) appears to be generally lower than that for non-inbred ones (F PED = 0).
To investigate this, we compared individuals from the two sets with virtually identical F values estimated by F ROH-0.1 (Table S7  arising from more ancient inbreeding, had likely been subjected to more intense purging than longer ROH belonging to F PED = 0.25 individuals, and arising from more recent inbreeding. This is in line with the observation in human data that long ROH are more enriched in deleterious mutations than short ones (Szpiech et al., 2013).
Simulation results obtained by Wang (2016) indicated that SNPby-SNP based measures offer accurate ID estimates as long as the number of markers is high enough (around 10,000). We did not find large differences when reducing the number of SNPs in the mean of SNP-by-SNP-based estimates of F and ID, although SEs increased considerably (Table S4, Supplemental File S3). However, ROH-based estimates were the most affected, as expected. Previous studies comparing the use of SNPs obtained from complete sequences or from SNP chips indicated that a large number of SNPs is required to properly detect short IBD segments (Ferenčaković et al., 2013;Purfield et al., 2012;Zhang et al., 2015). In agreement with Ferenčaković et al. (2013), we observed an increase of F ROH-1 for a decreasing number of SNPs up to 20,000, with a consequent decrease of ID. Whole genome sequencing provides information on rare alleles segregating at low frequencies, which can affect inbreeding estimates, particularly F YAN . Eynard et al. (2015) found significant differences on F YAN when obtained from common SNPs (MAF ≥0.05) or from rare alleles (MAF ≥0.01) using whole genome sequence data from Holstein bulls.
We did not find differences in F YAN estimates by applying a MAF pruning, probably because the low sample size considered, although we did observe a slight reduction in the estimate of ID.
As a concluding remark, our results support the use of molecular markers in inferring inbreeding depression, even in extreme situations with a minimal sample size of unrelated individuals derived from a single generation, and when the inbreeding measures are only partially associated with the measured traits.

ACK N OWLED G M ENTS
We thank two anonymous referees and the associate editor for useful comments on the manuscript. This work was funded by grants

CO N FLI C T O F I NTE R E S T S TATE M E NT
We declare no conflict of interests.