Prediction of fitness under different breeding designs in conservation programs

The expected change in fitness under inbreeding due to deleterious recessive alleles depends on the amount of inbreeding load harbored by a population, that is, the load of deleterious recessive mutations concealed in the heterozygous state, and the opposing effect of genetic purging to remove such a load. This change in fitness can be thus predicted if an estimate of the inbreeding load and the purging coefficient (the parameter that quantifies the amount of purging) are available. These two parameters can be estimated in pedigreed populations, as has been shown for populations under random mating. A question arises whether these parameters can also be estimated under other breeding systems as well as whether they allow accurate prediction of the corresponding expected change in fitness. In conservation programs, it is usually recommended to preserve genetic diversity by equalizing contributions from parents to progeny and avoiding inbred matings. Regular systems of inbreeding have also been proposed as breeding conservation strategies to purge the inbreeding load. Using computer simulations, we first test a method to jointly estimate the initial inbreeding load and the purging coefficient in populations subjected to equalization of parental contributions, circular mating, and partial full‐sib mating. Then, using the expected values of the inbreeding coefficient and the variance of family size, as well as the estimates of the inbreeding load and the purging coefficient, we make simple predictions of the change in fitness over generations under these breeding systems and compare them with simulation results. We discuss how these fitness predictions can help undertaking conservation designs under different breeding scenarios.


Introduction
Conservation programs are aimed not only at maintaining the highest possible levels of genetic diversity but also at avoiding inbreeding depression and population extinction (Frankham, Ballou & Briscoe, 2010). Inbreeding depression, the decline in fitness occurred because of the increase in inbreeding in a population, is a main determinant for the success of conservation programs (Frankham, 2005;Leroy, 2014;Hedrick & García-Dorado, 2016;Kariyat & Stephenson, 2019;Doekes, Bijma & Windig, 2021). Inbreeding depression is assumed to occur mainly by the expression in homozygosis of deleterious recessive alleles concealed in heterozygous state in noninbred individuals, the so-called inbreeding load (Crow, 1958;Roff, 2002;Charlesworth & Willis, 2009), although it may also be due to heterozygote advantage (Charlesworth & Willis, 2009;Charlesworth, 2015). When allele frequency changes due to selection are ignored, the rate of inbreeding depression equals the inbreeding load (Morton, Crow & Muller, 1956), which is generally expressed in terms of a number of lethal equivalents and varies widely depending on the species and particular populations (O'Grady et al., 2006;Plough, 2016).
The inbreeding load can be removed by genetic purging, the type of purifying selection that takes place when partial or fully recessive deleterious mutations are expressed in homozygosis by inbreeding (Gulisija & Crow, 2007;Hedrick & García-Dorado, 2016). The joint effect of inbreeding and purging on the change in population fitness can be described by an Inbreeding-Purging (IP) model developed by García-Dorado (2007. This model is based on a purged inbreeding coefficient, g, analogous to Wright's (1922) inbreeding coefficient (F) but corrected for the expected reduction in the frequency of recessive deleterious mutations due to purging. The purged inbreeding coefficient, in turn, depends on the purging coefficient, d, which measures the magnitude of purging and, for each locus, equals the recessive component of the deleterious effect that is concealed in the heterozygous state but is expressed in homozygosis. The efficiency of purging, that is, the proportional reduction of the frequency of deleterious alleles that purging is expected to cause in the long term (Pérez-Pereira et al., 2021a), depends on the rate of inbreeding (i.e., on the effective population size, N e ; Wright, 1931) and on the magnitude of d, purging becoming relevant for values of N e d > 1 ( García-Dorado, 2012).
If both the estimate of the inbreeding load harbored by a population and the purging coefficient ascribed to a particular breeding system are available, it is possible to predict the expected change in fitness over generations. The mean population fitness is expected to decline initially because of inbreeding depression but later recover partially or almost fully depending on the efficiency of purging. The joint estimation of both the inbreeding load and the purging coefficient and their use to obtain reliable predictions of the change in fitness over generations has been addressed for random mating (RM) populations (López-Cortegano et al., 2018). However, it is generally recommended that captive populations be managed to restrict inbreeding and maintain the highest possible genetic diversity Frankham et al., 2010). This is achieved by managing the choice of breeding adults to minimize their average coancestry (minimum coancestry contributions; Caballero & Toro, 2000) or, when the relatedness of the founder individuals is unknown, simply by equalizing parental contributions and avoiding mating between sibs (Caballero & Toro, 2000;Gulisija & Crow, 2007). Optimal contributions can also be applied to reduce the average coancestry while maximizing the genetic response in livestock breeding programs (Meuwissen, 1997;Meuwissen & Sonesson, 1998). One consequence of these management methods, where the contribution of offspring to the breeding pool is settled by the protocol, is that selection among families is depleted and there is, therefore, relaxation of the intensity of selection (Schoen, David & Bataillon, 1998;Fernández & Caballero, 2001). Thus, it is expected that the accumulation of inbreeding depression in managed captive populations will be slowed, but there will also be a restricted amount of genetic purging. In addition, if selection among families is avoided or nearly so, inbreeding depression and purging can be different for distinct components of fitness: selection for maternal fecundity is relaxed by equalization of parental contributions, whereas selection for viability is restricted to operate within families.
There has also been a discussion on the possibility of promoting purging in conservation programs, for example, by carrying out circular mating (CM) of half-sibs (Kimura & Crow, 1963;Theodorou & Couvet, 2010, 2015 or by forcing some proportion of mating between highly related individuals (Hedrick, 1994;Wang, 2000;Ávila, Amador & García-Dorado, 2010;Moreno et al., 2015). Purging by nonRM is generally more efficient than purging in small RM populations, although the efficiency of the former is diminished in populations of small size (Glémin, 2003).
Breeding systems promoting inbred mating would allow a higher genetic purging of the inbreeding load, but at the cost of higher initial inbreeding depression, with the corresponding increase in extinction risk (O'Grady et al., 2006;Caballero, Bravo & Wang, 2017;Cheptou, 2019). Thus, these methods could only be assumed for species with a high reproductive rate, being therefore not recommended for large mammals or birds. However, it can be considered for other animal species or plants able to cope more easily with the initial inbreeding effects.
The estimation of both the inbreeding load and the purging coefficient, which are the necessary parameters to predict the expected change in fitness of a population, can be obtained by the software PURGd (García-Dorado, Wang & López-Cortegano, 2016). This software requires data of pedigrees with measures of individual fitness, and its prediction accuracy has been tested assuming breeding individuals mated randomly and contributing a random number of breeding offspring to the next generation (López-Cortegano et al., 2018). It remains to be proved if appropriate estimations of the inbreeding load and the purging coefficient can also be obtained in the case of populations subjected to reproductive management, such as those where contributions from parents to offspring are equalized or where inbred matings are forced. With reliable estimates of these parameters, it would be possible to predict, using the IP model, the expected change in fitness under the above-mentioned breeding designs. This can be relevant from a practical point of view, as would imply that the possible outcome of a continuous breeding conservation strategy can be predicted. The purging coefficient and the consequent expected changes in fitness not only depend on the particular breeding systems carried out, but also depend on their degree of achievement. For example, equalization of contributions (EC) from parents to offspring can be fully achieved for a species with high reproductive performance, but probably not in full for a species with low reproductive capacity, where some couples can fail to produce the required number of offspring. Here, we carried out a computer simulation study with a double objective. First, to evaluate the reliability of the estimation of the inbreeding load harbored by a population and the purging coefficient with the software PURGd assuming different breeding systems of interest for conservation programs. Second, to evaluate if accurate predictions of the expected change in fitness over generations can be made for these breeding systems considering inbreeding depression and purging.

Materials and methods
The simulation procedure consisted of two steps, starting with the simulation of a large base population, then sampling a small number of individuals to establish small populations (lines hereafter) subject to different management systems. The output of these simulations was a pedigree file from each line including the fitness of individuals. The pedigree file was then used as input by the software PURGd to jointly estimate the inbreeding load (δ) and the purging coefficient (d), which were then used to predict the expected change in fitness under the IP model, together with the corresponding predictions under a model without purging (neutral model).

Simulation of the populations
First, a fully panmictic base population of N diploid individuals was simulated during 10 000 generations under the action of mutation, selection, and drift. Two population sizes differing in an order of magnitude (N = 1000 and 10 000) were assumed in order to produce two widely different inbreeding loads. The simulation parameters considered had been used by López-Cortegano et al. (2018) to test the predictive value of the IP model for populations with random contributions from parents and RM. Nonrecurrent deleterious mutation was assumed to occur at a rate λ = 0.2 per haploid genome and generation, a selection coefficient s was obtained from a gamma distribution with shape parameter β = 0.33 (mean s = 0.2), and a dominance coefficient h was obtained from a uniform distribution between 0 and e −ks (Caballero & Keightley, 1994), where k is a constant to obtain an average dominance coefficient value of h = 0.283. The fitnesses of any (biallelic) locus genotypes were 1, 1 − sh, and 1 − s for the wild-type homozygote, the heterozygote, and the homozygote for the deleterious mutant allele, respectively, with allele frequencies p and q. A set of 300 neutral loci were also assumed with the same mutation rate. All loci were assumed to be unlinked and the fitness of each individual was calculated multiplicatively among loci. The inbreeding load (δ) was obtained as the sum for all loci of s(1 − 2h)pq (Morton et al., 1956) which, if allele frequency changes by selection are ignored, equals the expected rate of inbreeding depression. After 10 000 generations, the base population was assumed to be under a mutation-selection-drift balance, and the number of segregating deleterious mutation loci was about 4000 (N = 1000) or 41 654 (N = 10 000). The observed inbreeding load obtained with these population sizes was δ = 2.52 and 6.23, respectively.
In a second step, a sample of N = 24 random individuals was taken from the base population to found a small line. This step was repeated to obtain a total of 10 000 randomly sampled lines per scenario considered (described below). These lines were maintained for t = N = 24 generations keeping the same mutation rate and mutational parameters as for the base population, but with some differences regarding the fitness model and the management of the matings and of the offspring contributed by parents to the next generation. Separate sexes and monogamous matings were always assumed in the lines unless indicated otherwise. For some simulations, new mutation was not considered in the lines in order to test the differences between simulations and predictions obtained by ignoring new mutation.
For the simulation of the lines, we considered two components of fitness, fecundity (number of offspring produced) and viability (probability of survival from birth to reproductive age) in order to use a more realistic model of fitness. We considered a model where half the fitness loci are assumed to control fecundity and the other half viability, and another where 1/3 of loci control fecundity and 2/3 control viability, following the estimates of inbreeding load for these two fitness components found by O'Grady et al. (2006) in mammals and birds. The individual fitness values were calculated separately for each component (W f and W v ), although the fecundity of males had no effect on the number of offspring, as fecundity was assumed to be a maternal fitness component. Overall fitness (W) was obtained as the product of both components. The number of offspring produced per family was calculated as 2KW f(mo) , where K is a constant representing the reproductive rate of the population and W f(mo) is the fecundity of the mother. Two contrasting values of K = 2.5 and 10 were assumed to evaluate typically low and relatively high reproductive rates. Each individual offspring survived with a probability equal to its viability, W v , with surviving individuals becoming part of the pool of potential breeders. The population size could vary among generations (N ≤ 24) if the total number of offspring available was not large enough, and the population could become occasionally extinct if N < 2.
Four management systems were considered:

Random mating and variable contributions
In each generation, parents were randomized and each pair formed contributed offspring to the next generation in proportion to the maternal fecundity. A maximum of 24 individuals were randomly sampled from the alive offspring to become the parents of the next generation and sex was then ascribed randomly (12 males and 12 females). Monogamous matings (each parent was only involved in a single mating pair) were always carried out in order to make a fair comparison with the other breeding systems explained below.

Equalization of contributions
Parents were randomized and pairs were formed avoiding full-sib matings. Each pair contributed two offspring to the next generation. If this was not possible because there were less than two alive offspring from the pair, the missing offspring was taken from another randomly taken pair with surplus progeny. Although the most efficient method to maintain neutral genetic diversity is to manage the choice of breeding adults to minimize their average coancestry (minimum coancestry contributions; Caballero & Toro, 2000), this procedure quickly converges to the EC protocol and, when the relatedness of individuals in the initial generation of management is unknown (as assumed here), minimum coancestry contributions reduces to equalization of parental contributions (Caballero & Toro, 2000).

Circular mating
Half-sibs were mated following a circular scheme so that individual i was mated with individual i + 1 to produce a male offspring, individual i + 1 with individual i + 2 to produce a female offspring, and so on, keeping this ordering in consecutive generations (see fig. 2 Crow, 1963). Thus, each individual contributed two (half-sib) individuals to the next generation. When there were no surviving offspring in the mating, the required offspring was obtained from the nearest mating pair in the circular scheme with a surplus of surviving offspring.

Partial full-sib mating
This was similar to RM but a given proportion of the matings was forced to occur between full sibs and the remainder mated at random. The final proportion of full-sib matings could be higher than this intended proportion due to chance, but also lower if the number of full sibs was not sufficient to cover the expected number of full-sib pairs. Contributions were again proportional to maternal fecundity. Two proportions of full-sib matings were tested, β = 0.2 and 0.5. As in the RM system, the surviving offspring from all pairs was randomized and a total maximum of 24 individuals were randomly sampled. Thus, the number of progeny from fullsib or nonfull-sib pairs was independently approximate Poisson distributed. Again, monogamous matings were carried out, and each parent was involved in just a single mating pair.

Estimation of the initial inbreeding load and the purging coefficient with PURGd
From each of the simulated lines, pedigrees were saved following the format required by the software PURGd (available at https://www.ucm.es/gfm/mecanismos-geneticos and https://gitlab.com/elcortegano/PURGd). These pedigrees included information (fecundity, W f , and viability, W v , expressed as the probability of survival) on all the individuals generated in each generation, living, and dead since otherwise inbreeding depression for viability would be underestimated. The initial inbreeding load (δ) and the purging coefficient (d), which are the main predictive parameters in the IP model (García-Dorado, 2012), were jointly estimated from each of these pedigrees following the numerical nonlinear regression method for untransformed fitness values included in PURGd . Estimates were first obtained with the default options and then with the option of estimating δ from individuals with null ancestral inbreeding, that is, F a = 0 (as defined by Ballou, 1997, F a is the average proportion of an individual's genome that has been exposed to inbreeding in at least one of its ancestors). The d and δ values obtained for each pedigree were then averaged over replicates. The value of d quantifies the amount of purging against deleterious recessive mutations, hence its name. For a single locus under selection, d equals the difference between the heterozygote and the mid-homozygote values for the fitness trait undergoing natural selection during the purging process. Thus, it is analogous to the definition of d for quantitative traits by Falconer & Mackay (1996, p. 109) but applied to fitness defined as the selection target during purging. Therefore, for a fully recessive lethal allele (s = 1, h = 0), the value of d is s(1/2 − h) = 0.5, and for a partially recessive deleterious allele with s = 0.1 and h = 0.2, its value is 0.03. However, if, for example, selection acts only on one sex, the d value would be halved. Similarly, d would be zero for fecundity under any management system ensuring that all parents contribute the same number of breeding offspring to the next generation, even if fecundity defined as number of offspring per parent is different for different individuals and undergoes inbreeding depression. Considering multiple loci with variable deleterious recessive effects, an effective purging coefficient should be used, which represents the d value that, when included in the IP model, produces the best fit for overall fitness (García-Dorado, 2012).
The mating systems studied differ mainly in the presence or absence of panmixia and in the contributions from parents to progeny (variable contributions for RM and PFS, and equal contributions for EC and CM). These differences could have an impact on the opportunities of purging of the inbreeding load, with possible different outcomes depending on the fitness component, fecundity, or viability.

Prediction of the change in fitness across generations
If estimates of the initial inbreeding load of the population and the purging coefficient are available, it would in principle be possible to predict the expected change in fitness over generations under any breeding system. The estimates of δ and d should ideally be estimated from the particular breeding system in question. However, it has particular interest to test whether the change in fitness can be accurately predicted by the model under different breeding systems if the estimates were obtained beforehand under the RM scenario. First, predictions of the expected N e and the expected inbreeding coefficients (F t ) of each generation (t) for the different breeding systems should be computed as explained in Appendix S1. Then, these predictions of N e and F t can be used, along with the available estimates of δ and d, to obtain the expected change in each fitness component separately and then to obtain overall fitness as the product of them. Thus, with the predictions of F t and N e for the different breeding systems, the purged inbreeding coefficient (g t ) for each generation t can be predicted as follows: (García-Dorado, 2012), where d is the purging coefficient expected under the pertinent breeding system, to be inferred as explained below from the estimate of d that was previously obtained under the RM scenario. Note that, as generations proceed, g t approaches an asymptotic value lower than one, where all the initial inbreeding load has been removed due to the combined effects of genetic drift and purging. If d = 0 (no purging), g t = F t . The larger is d, the smaller is the asymptotic g value and the role of genetic drift is reduced, with fewer deleterious alleles becoming fixed. To compute predictions for systems implying equalization of family contributions (EC and CM), we computed g t by using d = 0 in the case of fecundity and using the d estimate obtained under RM divided by √2 in the case of viability. This accounts for the fact that, in systems where the number of breeding offspring contributed per individual is settled by the management protocol, purging is expected to be relaxed for fecundity traits (i.e., d = 0) and is limited to operate upon the within lines additive variance for viability traits (García-Dorado, 2012) (the √2 correction appears because the selection effects are proportional to the root square of the variance, and only half of the additive variance is available to withinfamily selection). Then, the change in a fitness component (fecundity or viability) can be predicted as where W t is the expected mean fitness component at generation t, and W 0 and δ are the initial mean fitness component and inbreeding load, respectively, estimated also under the RM scenario. The change in overall fitness was predicted as the product of the predictions obtained for viability and fecundity. Under systems intended to equalize parental contributions, where this equalization was not fully achieved, this expression would predict too low fitness (too little purging).
We also obtained neutral predictions in order to describe the expected change in fitness in the absence of purging (d = 0 in equation 1). For this model, fitness predictions are obtained by W t ¼ W 0 Á e ÀδFt , where W 0 and δ correspond to the base population values observed in simulations.
In order to see the consequences of stopping the management process, for example, after the reintroduction of captive populations in the wild, the simulations were continued after the management period by assuming RM without management.

Results
Predictions of the variance of contributions from parents to offspring and the inbreeding coefficient Fig. 1 shows observed values of the variance of contributions from parents to progeny (S k 2 ) for each of the breeding systems and two reproductive rates (K). For CM and EC, the variance of contributions from parents to offspring should be S k 2 = 0 but, due to the failure of some families to produce two living descendants, the observed values were nonzero. For RM and PFS, S k 2 was also larger than expected, especially with a high reproductive rate (K = 10). Under a neutral model, however, values were very close to expectations for all breeding systems, as expected (Fig. S1).
Observed and predicted inbreeding coefficients are also shown in Fig. 1. As expected, the EC method produced lower inbreeding than RM, whereas methods promoting inbred matings (PFS and CM) produced higher F values. For the largest reproductive rate (K = 10), observed and expected values of F were generally close to each other for all breeding systems and there were almost no line extinctions (less than 0.1% of the replicates). For a low reproductive rate (K = 2.5), however, the values of F fell below predictions for breeding systems with forced inbreeding (PFS and CM), probably as a consequence of the higher extinction risk of lines under these systems (Fig. 2). Taking as reference the RM system, CM and PFS (particularly PFS.0.5) led to a dramatic increase of the extinction risk during the time scale of the simulations relative to that of RM, whereas EC led to some reduction. Thus, at the last generation, the extinction rates for CM and PFS.0.5 were 3.5 and 4.3 times larger than that for RM, respectively, and that for EC was 30% lower.

Estimates of the initial inbreeding load (δ)
The total inbreeding load (δ) obtained in the simulated base population of N = 1000 individuals was δ = 2.52 lethal equivalents (approximately half for each fitness component). The estimates of δ were generally uniform across mating systems (Fig. 3), although they were somewhat smaller than the true value. When the option of PURGd to estimate δ from individuals with F a = 0 was used, the average estimates were very close to the expected ones (Fig. 4).
Estimates of the purging coefficient (d) for each breeding system For RM and PFS, the estimates of d for fecundity were approximately half those estimated for viability systems because fecundity was simulated as a female trait, so purging would act only on half the population.
We will first consider results for scenarios of high reproductive rate (K = 10). In general, the values of d were not highly affected by the absence of panmixia as long as there were variable contributions from parents to breeding progeny, the d values being very similar for RM and PFS, although somewhat larger for the latter. On the contrary, strategies involving equalization of parental contributions (EC and CM) led to a relaxation of selection for fecundity giving lower values of the purging coefficient (Fig. 5a). Although theoretically, we would expect d to be 0 for fecundity in these cases, it was not so because equalization of parental contributions (S k 2 = 0) could not be fully achieved, as some families failed to produce progeny due to female sterility or low progeny viability (S k 2 > 0; Fig. 1).
Regarding viability, under complete equalization of parental contributions (S k 2 = 0¸d = 0 for fecundity), selection would act just upon the viability and only within families, thus using half the viability additive variance. Therefore, we could expect a reduction in d for viability in a proportion of 1/√2 with respect to RM. The estimated values (0.21 for CM and 0.25 for EC) were close to this theoretical expectation (0.23). When the reproductive rate was reduced to K = 2.5, estimates of d generally increased for all systems, but especially for CM (Fig. 5b). There are two possible explanations for this result, a greater purging both between and within families, or mainly purging between lines due to extinctions. In fact, extinctions were remarkable for CM and PFS.0.5 (Fig. 2). To evaluate the impact of such extinctions on the estimates of d, an alternative model of distribution of offspring was considered, where the total number of descendants was calculated per generation as 2KN. Such offspring was then distributed among families, with a probability of assignment proportional to the fecundity of the mother. Of the 1000-2500 simulated lines, no extinction was recorded under any management system or K. The estimates of d were similar between K = 10 and 2.5 (Fig. S2). Thus, most of the increase in the d values (and thus purging) for low reproductive rate was due to the extinction of lines.
The above results refer to the scenario in which fecundity was controlled by half of the loci and viability by the other half. When the corresponding proportions were changed to 1/3 and 2/3 (Figs S3 and S4), the results were not qualitatively different (cf. with Figs 3 and 5).

IP prediction of the change in fitness across generations
Predictions and observed values of fitness over generations for each breeding system are shown in Fig. 6. The predictions were made using the estimates of δ and d obtained from the RM scenario. As expected, predictions were close to observations for RM, but also for scenarios PFS.0.2 and EC (Fig. 6). The initial depression was predicted properly, especially with K = 10. For the systems with lower effective size (PFS.0.2 and PFS.0.5), there seemed to be a trend to overpredict long-term fitness (particularly under PFS.0.5), suggesting that part of the inbreeding load was due to alleles that have d < 1/N e in these lines. On the contrary, fitness was underpredicted under CM, severely for a low reproductive rate (K = 2.5). This latter occurred because the average fecundity was largely underpredicted for CM (Fig. S5), as the predictions assumed a null variance of contributions from parents to progeny (and, therefore, null purging) which could not be achieved in practice. Predictions for the viability component of fitness under CM were relatively accurate (Fig. S5). In fact, when the estimates of δ and d obtained from the CM breeding system were assumed to obtain these predictions, the predicted change in fitness was much closer to the simulated one (Fig. S5). Under a neutral model (only genetic drift, both during simulations and when computing predictions), there were also fitness overpredictions, which were partially corrected by avoiding new mutation in the simulations (Fig. S6), as the IP model does not consider mutation.
The above results assumed that the inbreeding depression in the base population is δ = 2.5. Results for a base population of size N = 10 000, for which the initial δ is 6.2, are shown in Appendix S2. Some quantitative differences were obviously found with respect to the previous ones, but the main conclusions could be reached. Thus, for example, the values of δ were well approximated, with about a maximum underestimation of 20% (Fig. A3 in Appendix S2). In addition, the estimated values of the purging coefficient (Fig. A4 in Appendix S2) were similar to those found in Fig. 5, although the values obtained with CM and EC were considerably larger as the larger inbreeding depression implied a poorer accomplishment of the management protocols even for K = 10 (see Fig. A1 in Appendix S2). Finally, the evolution of fitness was generally well predicted ( Fig. A5 in Appendix S2), except for a larger underestimation for the method CM.

Change in fitness after the management period
The different outcomes generated by the management with the different breeding systems are summarized in Fig. 7, Figure 4 Estimates of the inbreeding load (δ) of the base population referred to fecundity, viability, and overall fitness, obtained with the software PURGd with the option of estimating δ from individuals with F a = 0. Panels a and b indicate the results for different reproductive rates, K = 10 and 2.5, respectively. Boxes correspond to 50% of the estimates (interquartile range), and the mean is indicated by the horizontal line and the corresponding value indicated below. Whiskers indicate the minimum and maximum values within 1.5 times the interquartile range under the 25th percentile and over the 75th percentile, respectively. Horizontal dotted lines indicate the true value of δ obtained in simulations. CM, circular mating of half-sibs; EC, equalization of contributions avoiding mating between full sibs; PFS, partial fullsib mating at proportions 0.2 and 0.5; RM, random mating and variable contributions. which shows the overall fitness (W), inbreeding load (δ), and percentage of extinct lines at the end of the management process. Method CM showed a substantially lower fitness than the other methods for a reproductive rate K = 10, but there were no large differences in fitness for K = 2.5 between methods. Method CM, along with PFS.0.5, were those producing the highest extinction risk, whereas the lowest extinction risk but the highest inbreeding load was maintained by EC.
An interesting question is what happens if the management methods are stopped and RM without management is restored, for example, after the reintroduction of captive populations in the wild. The changes seen after five generations with the same population size and no management are shown in Fig. S7, and the overall changes in Fig. 7 by means of intervals and their corresponding arrows pointing out the direction and amount of the change. A large increase in overall fitness was observed for CM and a substantial drop in inbreeding load occurred for EC, but all methods suffered additional extinctions. In Fig. S7, we can observe that the large increase in fitness for PFS and especially CM occurred immediately after restoring panmixia.

Discussion
Inbreeding depression and purging of the deleterious genetic load are relevant and widespread phenomena with important consequences in the fields of evolution, conservation genetics, and animal and plant breeding (Lynch & Walsh, 1998, chap. 10;Caballero, 2020, chap. 8). Analyzing simulated Figure 5 Estimates of the purging coefficient (d) referred to fecundity, viability, and overall fitness and obtained with the software PURGd. Boxes correspond to 50% of the estimates (interquartile range), and the mean is indicated by the horizontal line and the corresponding value indicated below. Whiskers indicate the minimum and maximum values within 1.5 times the interquartile range under the 25th percentile and over the 75th percentile, respectively. Outliers are indicated by dots. Panels a and b indicate the results for different reproductive rates, K = 10 and 2.5, respectively. CM, circular mating of half-sibs; EC, equalization of contributions avoiding mating between full sibs; PFS, partial full-sib mating at proportions 0.2 and 0.5; RM, random mating and variable contributions. data with the IP model (García-Dorado, 2012), we show that the initial inbreeding load of a population and the purging coefficient ascribed to a given breeding system can be estimated rather accurately for a variety of designs that can either enhance or hinder genetic purging and that are recommended for the genetic management of endangered or productive populations. If an estimate of the initial inbreeding load (δ) and the purging coefficient (d) are available, we have shown that it is possible to predict the change in fitness over generations for these breeding systems considering the opposed effects of inbreeding depression and purging (García-Dorado, 2007. This can be particularly relevant for the design of conservation programs as it allows the possible outcome of the program to be foreseen. The expected change in fitness depends on the values of δ and d, and on the effective size (N e ) of the population as shown in fig. 1 of Pérez-Pereira, Caballero & García-Dorado (2021b) for RM populations, but it will also depend on the breeding system applied, as shown here. Thus, populations able to purge their inbreeding load will be less susceptible to the negative effects of inbreeding in the medium to long term, as long as they overcome the initial phase of inbreeding depression before the effects of purging become visible. Here we focused on the effects of the different management methods on fitness, as this was our main objective. Genetic purging, as background selection in general, implies a reduction of genetic variation (Charlesworth & Charlesworth, 2010). However, it is convenient to note that purging may slow the loss of genetic diversity compared with that expected in the absence of purging when there is tight linkage, as shown theoretically (Wang et al., 1999;Bersabé et al., 2016). The estimates of δ obtained by PURGd were close to the true simulated value but showed a systemically downward bias (Fig. 3). This was expected to some extent based on previous results. Using simulated data, López-Cortegano et al. (2018) showed that longer periods of slow inbreeding give higher power to detect purging, which is then more efficient according to the theory, but is also accompanied by a greater underestimation of both δ and d when they are estimated jointly. Our pedigrees, encompassing t = N generations, were long enough to detect purging while inducing a relatively small underestimate of δ. This agrees with similar purging estimates obtained in ungulate species, where purging was detected in pedigrees with about t = 9 generations with overall effective population sizes N e slightly over 10, while detection of purging failed when the number of generations (t = 7) was well below the effective size (N e = 39) (López-Cortegano, Moreno & García-Dorado, 2021). Furthermore, an effective size below 10 can also be too small to allow for purging efficient enough to be detected. This can explain why many studies of real data have failed to detect purging as, for example, the meta-analysis by Boakes, Wang & Amos (2007), or the recent analysis of Kangaroo rat data using PURGd, as discussed below (Willoughby et al., 2019). When δ was estimated using only individuals with null ancestral inbreeding and, therefore, no purging (F a = 0), the mean estimated δ values were almost identical to the true simulated one (Fig. 4). Thus, if the initial δ value is to be estimated with no bias, this should be done considering only individuals with null ancestral inbreeding (F a = 0).
There are other alternative methods to PURGd to estimate inbreeding depression depending on the type of data available. Nietlisbach et al. (2019) compared different statistical models to estimate the inbreeding load, showing that maximum likelihood methods (Kalinowski & Hedrick, 1998) produce unbiased estimates. Willoughby et al. (2019) applied both the Kalinowski & Hedrick (1998) method and the IP model with PURGd to estimate the inbreeding load of a wild population of banner-tailed kangaroo rats using a reconstructed pedigree. Estimates of the inbreeding load ranged from 0.98 to 4.7 haploid lethal equivalents depending on the fitness trait, and these were very similar for the two estimation methods. Inbreeding depression can also be estimated in the absence of pedigree records using molecular markers (Howard et al., 2017). In fact, molecular markers can be more powerful to detect inbreeding depression than pedigree records if they are numerous (more than 10 000; Wang, 2016;Nietlisbach et al., 2019). However, the estimation of inbreeding depression from molecular markers has also difficulties and estimates can be biased and noisy (see, e.g., Yengo et al., 2017;Nietlisbach et al., 2019;Caballero, Villanueva & Druet, 2021). Nevertheless, it is not an objective of the present study to compare the estimation properties of different estimators of the inbreeding load, but to note that if these estimates and those of the purging coefficient for a given breeding system can be obtained, the change in fitness over generations can be predicted. Some alternative methods have also been developed to quantify the amount of genetic purging, either using regression models based on an ancestral inbreeding coefficient, F a (Ballou, 1997;Boakes et al., 2007) or defining several indices based on the opportunity of purging using pedigree information (Gulisija & Crow, 2007, as estimated in the reimplemented R package of PURGd -purgeRby López-Cortegano, 2021). However, only the IP approach allows an estimate of the predictive purging parameter d to be obtained and used to predict the consequences of purging on fitness and inbreeding load. Our results show that the d estimates are rather variable across simulated lines for any given breeding system (Fig. 5), pointing toward high stochasticity attached to this parameter. This variability can be mainly attributed to the different impacts of purging and drift among replicates, as a primary consequence of the random sampling of a small number of individuals from a large population where recessive deleterious alleles are expected to segregate at low frequencies. Part of the variability could also be ascribed to estimation error, but this is likely to be minor. For breeding systems equalizing parental contributions (CM and EC), we would expect purging to be absent for fecundity under a neutral model, which would lead to d = 0. Likewise, for viability, we would expect a reduction of the family component of the additive variance for viability. This led to estimates of d that were above those expected under perfect equalization (above 0 for fecundity and above the 1/√2 times the d estimates obtained under RM for viability). Even so, the values of d were substantially lower than that under RM. This result is in line with previous studies showing low statistical power to detect inbreeding depression when inbreeding is avoided to maximize genetic diversity (Kalinowski & Hedrick, 1999). In this scenario, both inbreeding depression and purging are expected to be reduced.
Several studies have obtained estimates of the purging coefficient from RM populations. Bersabé & García-Dorado (2013) obtained values of the overall purging coefficient against nonlethal deleterious alleles of about d ≈ 0.02 for small populations of Drosophila melanogaster with size N e < 12, whereas López-Cortegano et al. (2016) obtained values of about d ≈ 0.3 (or 0.2 for nonlethal mutations) under more competitive conditions for large populations (N e > 1000) of D. melanogaster, which are comparable to those simulated here. Willoughby et al. (2019) obtained nonsignificant estimates of the purging coefficient using PURGd, ranging from 0.006 to 0.014, which could be due to low experimental power related to the limited depth of the pedigree. Their reported average inbreeding coefficient was F = 0.07 and the reconstructed pedigree had a depth of six generations, suggesting an effective population size of around 35, so that six generations can represent a too short period for efficient purging effects to become appreciable. In addition, purging could have occurred previously, as noted by the authors. Another explanation can rely on the particular past demography of the population. If this was maintained under similar conditions and similar population size over a period of time much longer than the generations evaluated, the population may have reached a mutation-selection-drift equilibrium, where δ is constant and the fitness decline can be the result of the fixation of deleterious alleles and not inbreeding depression. Finally, as noted above, López-Cortegano et al. (2021) recently estimated large and significant d values in two captive gazelle populations that were maintained for about eight generations with N e > 10, and it is worthwhile to note that, in this latter case, management was intended to minimize coancestry between mating individuals but the contribution of offspring per female was not under control. Therefore, the available evidence suggests that purging can substantially reduce the fitness inbreeding depression in not too small populations (say, N e > 10) as far as the managing protocols allow for natural selection, as shown in our simulation study. In fact, recent experimental work in D. melanogaster shows that in large RM populations (N e > 1000), virtually all the initial inbreeding load can be purged if enough time is allowed (Pérez-Pereira et al., 2021a), which suggests that purging under RM is effective in large populations.
As mentioned above, a prediction of the change in fitness considering the consequences of purging can be useful to predict the future outcome of a given breeding design in conservation programs as well as the purging expectations in comparison to neutral conditions. Using the estimates of δ and d obtained from the RM scenario, we obtained predictions for the change in fitness for the different breeding systems with a generally good fit (Fig. 6). This suggests that the values of δ and d estimated with RM are good enough to make predictions for the other breeding designs if account is taken of the reduced purging expected for the designs where parental contributions are equalized (EC and CM). It also shows that predictive equations to obtain the variance of family size and the expected coefficient of inbreeding for the different designs (Appendix S1) are rather accurate. accounting for mutation and nonpurging selection (García-Dorado, 2012), can be satisfactorily used in these conditions. In the case of CM, where each individual should contribute two offspring, the fitness average increased well above the IP prediction. This should be ascribed to perfect equalization not being achieved and, in the K = 2.5 case, to betweenlines selection associated with the large extinction risk experienced by these lines, where the initial rate of inbreeding is very high. Under partial full-sib (PFS) mating, the inbreeding depression was lower than expected, but purging seemed also to be either lower or less efficient. However, it is worth noting that, in these cases, the observed fitness of neutral lines was initially higher than expected and became lower later on. This suggests that the reason for the poor fit to predictions under PFS is that the increase of inbreeding was not accurately predicted.
For low reproductive potential, short-term extinction risk can be large and should be a main concern. Therefore, the methods generating higher levels of inbreeding (CM and PFS.0.5) produce the highest extinction risk, whereas the lowest extinction risk was maintained by EC at the cost of a higher inbreeding load (Fig. 7). Thus, as already noted before (Schoen et al., 1998;Fernández & Caballero, 2001), EC delays the increase in inbreeding but relaxes natural selection, so that less inbreeding load is purged and it is likely that extinction would become larger than under the other systems in the long term or after reintroduction (see below). In contrast, methods that force inbred matings produce more purging, but at the cost of higher inbreeding depression and extinction rate. It has been shown theoretically that forcing inbreeding by nonRM can induce more efficient purging than inbreeding induced by the small population size alone (Glémin, 2003). In accordance with this prediction, we observed a lower final inbreeding load for the nonrandom PFS mating schemes than for RM, but again with a higher extinction risk under low reproductive potential. Thus, from a practical point of view, a warning should be noted here, as not all populations will be able to deal with increased inbreeding depression unless the reproductive rate, or population size, is high enough, as previously noted (Theodorou & Couvet, 2015;Caballero et al., 2017). For instance, reproductive rates in mammals and birds are commonly below 10 (measured in terms of total descendants per year; Vance, Fahrig & Flather, 2003;Rytwinski & Fahrig, 2011;Quesnelle, Lindsay & Fahrig, 2014), and values below 3 can be found for which extinction due to increased inbreeding depression by a particular management method would be the main concern.
After a period of no management, for example, after reintroduction of the captive managed population in the wild, a large increase in fitness was observed in our simulations for CM due to restauration of Hardy-Weinberg frequencies and a substantial decline in inbreeding load occurred for EC (Fig. 7). Thus, forcing or avoiding inbred matings and equalizing or not contributions from parents will depend on the particular situation regarding the reproductive rate of the species and its capacity to circumvent the negative effects of inbreeding. For example, for large reproductive rates, where extinction is not an immediate concern, PFS and CM can lead to high average fitness after reintroduction with minimum inbreeding load, favoring population resilience to future demographic challenges, whereas equal contribution systems would lead to the opposite results. On the contrary, under low reproductive potential, PFS and CM are always to be avoided, as they are associated with too high immediate extinction rates, while both RM and EC lead to low extinction and similar fitness after reintroduction, although the latter gives larger inbreeding load which implies lower resilience to future inbreeding.
The simulation model used herein assumed mutational parameters for deleterious mutations within the ranges obtained empirically from mutation-accumulation experiments (Caballero, 2020, p. 156). We disregarded the introduction of favorable mutations, as this would be expected to have a negligible contribution in short periods of time and small populations, as assumed here. With the mutational parameters considered, and assuming a base population with size N = 1000 individuals, the expected inbreeding load of the base population reached a value (δ = 2.5) within the range of those usually observed in captive populations (Ralls, Ballou & Templeton, 1988) and in experimental populations of D. melanogaster (Ávila et al., 2010;López-Cortegano et al., 2016). If a ten times larger base population was assumed (N = 10 000 individuals), the observed inbreeding load amounted to δ = 6.2, a value observed for fitness in mammals and birds by O'Grady et al. (2006). However, the results obtained under this last scenario (Appendix S2) did not differ qualitatively from the former, and the same overall conclusions were obtained.
With RM and variable contributions (RM), the estimated purging coefficient considered was about d = 0.2, a value within the range of those obtained experimentally by López-Cortegano et al. (2016) in Drosophila populations for nonlethal alleles, as mentioned above. This value is taken as a reference to be compared with the amount of purging associated with other mating systems.
We assumed in our simulations and predictions that inbreeding depression is only due to the expression of (partially) recessive deleterious mutations. This is the scenario for which numerous estimates of mutational parameters are available (e.g., Caballero, 2020, p. 156) and thus its contribution to inbreeding depression can be modeled. It is possible that heterozygote advantage for fitness also contributes to inbreeding depression, but its relative contribution is at the moment fully speculative. According to most evidence, the contribution of overdominance to inbreeding depression is probably minor Roff, 2002;Charlesworth & Willis, 2009;Hedrick, 2012;Yang et al., 2017). In their genomic study with maize lines, Yang et al. (2017) did not neglect overdominance as a possible contributor to heterosis on some occasions, but they argued that a model of deleterious recessive alleles may largely explain not only the results for fitness-related traits but also for hybrid vigor. Some recent theoretical analyses of data in Drosophila, however, suggest that balancing selection could contribute to genetic variation in this species  (2015) was based on some arguable assumptions, such as mutation-selection balance, and other explanations rather than overdominance for fitness can be considered, while the results of Sharp & Agrawal (2018) could be affected by the bias of predictions ascribed to natural selection. Ayroles et al. (2009) analyzed the mode of action of genes differentially expressed in highly depressed inbred lines of D. melanogaster. They found that about 75% of genes presumably involved in inbreeding depression were additive or partially dominant, but 25% showed a pattern of overdominance. Similar studies found, however, much lower proportions of genes showing patterns of overdominance (<5%) (Gibson et al., 2004;Hughes et al., 2006). Even so, as acknowledged by Ayroles et al. (2009), the proportion of putative overdominant loci could also be ascribed to pseudooverdominance generated by linked deleterious mutations (see Schou et al., 2017;Waller, 2021), a likely explanation considering the reduced genetic length of the Drosophila genome and its multiple inversions. In addition, the trait analyzed by Ayroles et al. (2009) was the amount of gene expression, and it has been shown that there is not a clear relationship between differential gene expression due to inbreeding and inbreeding depression (Schou et al., 2018), so the mode of gene action analyzed by Ayroles et al. (2009) might not be related at all with inbreeding depression for fitness. In a large meta-analysis of estimates of selection coefficients in natural populations carried out by Thurman & Barrett (2016), 3416 estimates were ascribed to directional selection and only 70 to heterozygote advantage. Estimates can be very noisy and the small proportion of overdominant effects could be false positives, so that it is very difficult to extract clear-cut conclusions. Finally, it should be noted that the inbreeding load generated by heterozygote advantage for fitness cannot be purged. A recent study showed that the inbreeding load of two large populations of Drosophila (N e > 1000) was fully exhausted by genetic purging (Pérez-Pereira et al., 2021a), a result which cannot be accounted by genetic drift alone given the large effective population sizes. This removal of the inbreeding load was shown to be parsimoniously explained by genetic purging and drift operating on genetic variation for fitness due to partially recessive deleterious mutations (Pérez-Pereira et al., 2021a). In summary, although there might be a contribution from fitness overdominance to inbreeding load, this cannot be reliably quantified. Thus, the most conservative model to be assumed, as we did in this study, is that inbreeding depression arises mainly from partial dominance.
Our simulation results assume for simplicity that the environment is stable over generations. This is a reasonable assumption during a captive breeding period, as implied in most of our results. If the managed population is reintroduced into the wild, it is expected that the fitness effects of past purging will change, as the expression of the inbreeding load depends on the environment (Kristensen, Loeschcke & Hoffmann, 2008;Cheptou & Donohue, 2011;Reed et al., 2012;Pemberton et al., 2017). For example, no purging has been detected for specific traits that are not expressed during the purging process as heat shock resistance (Bijlsma, Bundgaard & Van Putten, 1999;Bundgaard et al., 2021). However, even though adaptation to captive conditions is likely to imply some detriment to adaptation in the wild, this trade-off can depend on just a small fraction of the mutations affecting fitness which do not need to make substantial contributions to the fitness inbreeding load. In fact, most deleterious mutations are expected to have some deleterious effect both in captive and wild conditions, these effects being usually larger in harsh environments (Halligan & Keightley, 2009). In agreement with this, the inbreeding load is usually larger when assayed in the wild (Keller & Waller, 2002;Armbruster & Reed, 2005;Fox & Reed, 2011), and purging under competitive conditions has been found to reduce the inbreeding depression expressed in noncompetitive ones . Therefore, although more results are needed in this respect, purging that occurred during ex situ conservation programs should in principle be expected to reduce the inbreeding load and inbreeding depression expressed after reintroduction into the wild (Swindell & Bouzat, 2006). In the final five-generation period without management considered in our simulations, we assumed, however, that the efficiency of purging is the same as for the management period. The reason is that our intention was to show the consequences of stopping the management breeding without other confounding factors.
Summarizing, the IP model offers a tool able to jointly estimate both the initial inbreeding load harbored by a population and the corresponding purging coefficient expressed under different breeding systems, involving the management of offspring contribution and/or mating and implying different intensities of purging. These parameters, along with calculations of the expected inbreeding coefficient for each breeding system, allow the prediction of the change in fitness over generations considering inbreeding depression and purging, which may help to foresee the possible outcomes of alternative conservation programs.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1. Left-hand graphs: Observed (columns) and expected (black thick lines) variance of contributions from parents to breeding progeny (S k 2 ). Right-hand graphs: Inbreeding coefficient (F) observed in simulations (dotted lines) and predictions (solid lines) for each mating system. Figure S2. Estimates of the purging coefficient (d) referred to fecundity, viability, and overall fitness, obtained with the software PURGd for simulations with an alternative model of distribution of contributions (the total number of offspring per generation is calculated as 2KN and is distributed among families with a probability of assignment proportional to maternal fecundity). Figure S3. Estimates of the inbreeding load (δ) of the base population referred to fecundity, viability, and overall fitness, obtained with the software PURGd. Figure S4. Estimates of the purging coefficient (d) referred to fecundity, viability, and overall fitness and obtained with the software PURGd. Figure S5. Change in fecundity (W f ), viability (W v ), or overall fitness (W) in the case of circular mating for different reproductive rates, K = 10 and 2.5. Figure S6. Change in overall fitness (W) in the absence of new mutation. Figure S7. Change in overall fitness (W) during 24 generations of management, followed by five generations of random mating without management (RM), for different reproductive rates (K = 10 and 2.5).
Appendix S1. Prediction of the expected inbreeding coefficient (F) and effective population size (N e ) for different mating designs.
Appendix S2. Simulation results for a base population with size N = 10 000 individuals and an initial inbreeding load of δ = 6.23.