Heterosis, Marker Mutational Processes and Population Inbreeding History

Corresponding author: Anne Tsitrone, CEFE-CNRS, 1919 Rte. de Mende, 34293 Montpellier Cedex 5, France., tsitrone{at}cefe.cnrs-mop.fr (E-mail)

Communicating editor: D. CHARLESWORTH

	ABSTRACT

TOP ABSTRACT THEORY RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C LITERATURE CITED

Genotype-fitness correlations (GFC) have previously been studied using allozyme markers and have often focused on short-term processes such as recent inbreeding. Thus, models of GFC usually neglect marker mutation and only use heterozygosity as a genotypic index. Recently, GFC have also been reported (i) with DNA markers such as microsatellites, characterized by high mutation rates and specific mutational processes and (ii) using new individual genotypic indices assumed to be more precise than heterozygosity. The aim of this article is to evaluate the theoretical impact of marker mutation on GFC. We model GFC due to short-term processes generated by the current breeding system (partial selfing) and to long-term processes generated by past population history (hybridization). Various mutation rates and mutation models corresponding to different kinds of molecular markers are considered. Heterozygosity is compared to other genotypic indices designed for specific marker types. Highly mutable markers (such as microsatellites) are particularly suitable for the detection of GFC that evolve in relation to short-term processes, whereas GFC due to long-term processes are best observed with intermediate mutation rates. Irrespective of the marker type and population scenario, heterozygosity usually provides higher correlations than other genotypic indices under most biologically plausible conditions.

THE existence of correlations between individual genotypes at marker loci and fitness-related traits has caused much debate among evolutionary biologists. Such correlations were initially used as an argument in favor of selection acting on the maintenance of allozyme polymorphisms in the controversy that has historically opposed selectionists and neutralists (DAVID 1998 ). Allozymes have been used for decades to detect correlations between multilocus heterozygosity (the number of heterozygous marker loci per individual) and fitness-related traits such as growth, viability, or physiological parameters (reviewed in MITTON and GRANT 1984 ; DAVID 1998 ). Such positive heterozygosity-fitness correlations (HFC) have been reported for various organisms, including marine bivalves (ZOUROS et al. 1988 ), salmonid fishes (LEARY et al. 1983 ), and pine trees (LEDIG 1986 ). HFC have recently been reported using restriction fragment length polymorphism (RFLP) markers (POGSON and FEVOLDEN 1998 ) and microsatellite markers (BIERNE et al. 1998 , BIERNE et al. 2000A ).

The observation of significant HFC with noncoding DNA markers makes it clear that at least some of the correlations are not due to direct effects of the marker genes on the phenotype. Associative overdominance refers to any kind of HFC not due to a direct effect of marker genes, but to genetic associations between the markers and fitness genes (DAVID 1998 ). The first kind of association is linkage desequilibrium due to genetic drift (correlation of allelic state within gametes). This has been identified as a possible cause of HFC in small populations when there is physical linkage between fitness genes and marker genes (OHTA and KIMURA 1970 ; PAMILO and PALSSON 1998 ). The second kind of association, identity disequilibrium due to variance in inbreeding (correlation of homozygosity between loci across the whole genome), has been identified as a major source of HFC in several theoretical models (OHTA and COCKERHAM 1974 ; CHARLESWORTH 1991 ; ZOUROS 1993 ). Here we focus on this second process. Variance in inbreeding generates HFC because more inbred individuals are both more homozygous for their marker loci and less fit due to inbreeding depression. HFC may thus be a powerful tool to analyze the fitness consequences of inbreeding (CHARLESWORTH and CHARLESWORTH 1999 ; PEMBERTON et al. 1999 ). However, inbreeding itself may be caused by very different population processes. The most obvious is the mating system, which generates "short-term" inbreeding, i.e., inbreeding caused by one or a few generations of consanguineous matings. This could explain HFC in large, partially selfing populations of pine trees (LEDIG 1986 ; BUSH and SMOUSE 1991 ). "Long-term" inbreeding, on the other hand, involves both recent coalescence events and coalescence events deeper in the pedigree. For example, when two isolated populations come into contact, hybrid offspring are more "outbred" than nonhybrid offspring. In this case, inbreeding may have built up during a long history of isolation of the parental source populations. Such a scenario was invoked to explain the HFC detected in the red deer population of the Isle of Rum (COULSON et al. 1998 ) and in the harbor seal population breeding on Sable Island (COLTMAN et al. 1998 ; PEMBERTON et al. 1999 ).

Previous models of associative overdominance implicitly neglect marker mutation, because they focus on short-term inbreeding and on markers with low mutation rates (e.g, allozymes). Now that highly mutable markers (i.e., microsatellites) and long-term scenarios are included in HFC studies, a theoretical assessment of the importance of mutation is needed. This approach was implicitly followed by COULSON et al. 1998 who proposed the use of a new individual genotypic index, rather than heterozygosity, to account for the marker mutation properties (in their case, microsatellites). The model usually assumed for microsatellites is stepwise mutation (one repeat unit added or removed in each mutation event; OHTA and KIMURA 1973 ; VALDES et al. 1993 ). This mutation model suggested the definition of an index, d², the squared difference in repeat units between the two alleles of an individual (COULSON et al. 1998 ), whose distribution is closely related to the distribution of coalescence times under such a mutation model (PRITCHARD and FELDMAN 1996 ). Empirical studies have sometimes found d² to correlate with fitness traits in samples where heterozygosity does not correlate significantly with these traits (COLTMAN et al. 1998 ; COULSON et al. 1998 ). It has thus been suggested that heterozygosity is suitable for detecting short-term inbreeding, whereas d² provides additional information about long-term inbreeding, due to the mixture of formerly isolated subpopulations (PEMBERTON et al. 1999 ; COULSON et al. 1999 ; MARSHALL and SPALTON 2000 ).

In this article we assess the significance of these arguments using a theoretical approach. The influence of marker mutation on genotype-fitness correlations (GFC) due to inbreeding is investigated by comparing (i) different mutation models and mutation rates, (ii) different population inbreeding histories involving various time scales, and (iii) different genotypic indices related to the mutational processes of the marker.

	THEORY

TOP ABSTRACT THEORY RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C LITERATURE CITED

The model:
The rationale of the model is as follows: Population history can generate variance in inbreeding among individuals, depending on individual pedigrees. The population can thus be partitioned into "inbreeding classes." First, inbreeding levels are associated with genotypes at neutral markers. The latter can be obtained from identity-in-state (IIS) relationships among marker alleles for a given marker mutation model and population scenario. Second, in the presence of inbreeding depression, the inbreeding level determines the value of the fitness trait. Thus, individual phenotype and genotype correlate through individual variation in the inbreeding level. For each mutation model, genotypic index, and population scenario assumed, the correlation coefficient (X, W) between X (a given index) and the fitness trait W is derived. As

(1)

three moments must be computed: the covariance of the fitness trait and X, cov(W, X), and the variances ²(W) and ²(X). Below, we describe how analytical expressions for these moments can be obtained. These computations rely on IIS probabilities that depend on the mutation model assumed. Mutation models for the marker genes and corresponding indices are described in the next section.

Mutation models and individual genotypic indices:
In a second step we develop theoretical models that approximate mutational processes at DNA markers, such as microsatellites, RFLP markers, or neutral DNA sequences. Only the first two categories of markers have been used in GFC studies, but the third kind might also be used in the future.

Stepwise mutation model: Under a strict stepwise mutation model (SMM; OHTA and KIMURA 1973 ) with mutation rate u, an allele with i repeat units is assumed to mutate only to the i - 1 or the i + 1 states, each with probability u/2 per generation. This model is classically taken as an approximation for mutation at microsatellite marker loci, for which u ranges from 10^-6 to 10^-2 per generation (JARNE and LAGODA 1996 ; ESTOUP and ANGERS 1998 ). Two indices are considered. The first is individual heterozygosity H, which takes the value 0 when the marker locus is homozygous and 1 when it is heterozygous. The second is d², the squared difference in repeat units between the two alleles of an individual (COULSON et al. 1998 ).

K-alleles model: Under a model first formulated by CROW and KIMURA 1970 , there is a finite number K of possible allelic states, and each allele can mutate to any other at rate u/(K - 1). A two-allele model can be taken as an approximation for mutation at an RFLP locus, as only two alleles ("cut" and "uncut") need to be distinguished, and for single nucleotide polymorphisms (SNPs; KUHNER et al. 2000 ). We assume the mutation rate of these DNA markers to be fairly low, even though no direct data are available. Estimates of the substitution rate per nucleotide site and per generation are typically of the order of 10^-8 or less, depending on the organism under consideration (LI 1997 ; DRAKE et al. 1998 ; NACHMAN and CROWELL 2000 ). For a SNP site, u should thus be at most 10^-8 per generation. Assuming that an RFLP polymorphism is due to nucleotide substitutions rather than to indels and that the RFLP locus is composed of few nucleotides, the mutation rate for an RFLP should be of the order of 10^-8 or 10^-7 per locus per generation. Although mutational processes affecting RFLP loci may be asymmetric, destroying any particular polymorphism more often than recreating it, we do not expect our model to be very sensitive to this asymmetry, given that mutation rates will be low (but this remains to be tested). The only index studied for this kind of marker locus is individual heterozygosity.

Infinite-alleles and infinite-sites models: Under the infinite-alleles model (IAM; KIMURA and CROW 1964 ), each mutation event produces a new allele at rate u per generation. The infinite-sites model (ISM; KIMURA 1969 ) shares the same property but is more specifically designed for DNA sequences. In this model the number of nucleotide sites in the sequence is assumed to be so large that each new substitution occurs at a site that has not mutated before. The total mutation rate of the sequence is u = lµ, where l is the sequence length in base pairs and µ is the substitution rate per nucleotide site per generation (see above). The first index considered is sequence heterozygosity H, which takes the value 0 when the two alleles of the marker are strictly identical in sequence and 1 when there is at least one different nucleotide site between the two allelic states. Note that analytical derivations of heterozygosity under the ISM and the IAM with mutation rate u are strictly equivalent. For neutral DNA sequences, the number of nucleotide differences between the two alleles of an individual, denoted by p, is also used as an alternative index.

The partial selfing model and single-locus heterozygosity-fitness correlation:
We investigate a simple scenario in which HFC relates to the actual mating system of the population. This scenario examplifies short-term inbreeding. Consider a large population of size n at inbreeding equilibrium, with freely recombining loci and one marker locus with a given mutation model. A proportion S of offspring is produced by selfing at each generation, whereas a proportion 1 - S comes from outcrossing events. Each individual is characterized by the number of generations of selfing in its pedigree, k, starting from the most recent outcrossing event (0 k ). The population is thus partitioned into inbreeding classes C_k, 0 k , each consisting of individuals having the same inbreeding level (i.e., the same k). Note that, under this model, there is an infinite number of inbreeding classes, although in practice only classes with low k are largely represented. In the absence of selection, inbreeding class k has frequency Pr(C_k) = (1 - S)S^k. Selection against homozygous genomes reduces the frequency of classes with high k. The effect of selection thus resembles a reduction in the selfing rate S to a lower value S_sel that can be computed as explained in DAVID 1999 (Appendix A). In practice, DAVID 1999 shows that frequencies of inbreeding classes computed using S_sel instead of S provide a good approximation of the frequencies with selection. In what follows, we therefore take into account selection against homozygous genomes simply by replacing the "raw" S value by the value of S_sel.

Let us define the probabilities of IIS within a population: Q₀ for a pair of genes drawn from the same individual and Q₁ for genes from different individuals. In what follows, moments are expressed as functions of Q₀ and Q₁ that are derived in Appendix A under each mutation model. The inbreeding coefficient in class k, f_k, is defined by

(2)

where Q_|k is the probability of identity of genes in an individual of class C_k. Under a general inbreeding load model assuming additive effects of fitness genes (MORTON et al. 1956 ) the value of a fitness trait W in a C_k individual can be written as

(3)

where W₀ is the value of the fitness trait for outbred individuals (i.e., belonging to C₀); ß is the inbreeding load (i.e., the fitness reduction among completely inbred individuals as compared to outbred individuals); and is a random effect with mean 0 and variance ², assumed to be independent of k. Moments of the fitness trait W are, based on Equation 3,

(4a)

(4b)

where E(f) and ²(f) are the mean and variance of f_k over all k-values.

Moments of heterozygosity H are simply expressed as functions of probabilities of IIS:

(5a)

(5b)

Within a given inbreeding class C_k, heterozygosity is not correlated with the fitness trait so that, conditioning on the inbreeding class k,

(6a)

where E(H|C_k) = 1 - Q_|k, and E(W|C_k) = W₀ - ßf_k. First- and second-order moments of f_k and Q_|k are needed to compute Equation 4aEquation 4b. Neglecting mutation since the most recent outcrossing event, 1 - Q_|k , so that f_k = 1 - ()^k. This yields

and

where F_IS is the classical F-statistic (WRIGHT 1951 ). Note that F_IS can be estimated directly using the heterozygote deficiency (WEIR and COCKERHAM 1984 ). Finally,

(6b)

Note that cov(H, W) is null for S = 0 and S = 1, since there is then no variance in the inbreeding level among individuals. Equation 4b, Equation 5b, and Equation 6b can now be used to derive the correlation coefficient between heterozygosity and the fitness trait from (1). The maximum correlation coefficient (_max), assuming no within-inbreeding class variance for the fitness trait (² = 0 in Equation 3), is

(7a)

²_max does not depend on the inbreeding load ß, but only on the marker genetic diversity 1 - Q₁ (see Appendix A for an analytical expression of Q₁ under each mutation model) and on the selfing rate S. In the presence of within-pedigree variance for the fitness trait (² > 0), ² is

(7b)

Extension of the partial selfing model to other genotypic indices:
We first focus on the squared difference in repeat units for a microsatellite marker under the stepwise mutation model of mutation. Cov(d², W) is derived as in (6a), replacing H by d². Moments of d² are obtained as detailed in Appendix A (Equation A7). This yields

(8) and thus

(9) and ² is as in (7b).

We then focus on another genotypic index: the number of nucleotide differences for a neutral sequence under the infinite-sites model of mutation. The correlation coefficient between p and W under the ISM is computed similarly to the correlation coefficient between d² and W. Moments of p are derived in Appendix A (Equation A5) and Cov(p, W) is as in (8). We obtain

(10)

and ² is as in (7b).

The population admixture model (hybridization):
We next consider a simple situation where the HFC is related to long-term processes. We assume a large, random mating population of size N at mutation-selection equilibrium, which splits into two randomly mating finite subpopulations of size n << N. Let the two subpopulations diverge for a long time (assumed to be equal to N generations for the sake of simplicity) without any gene flow between them and then merge into a single, infinite, panmictic population. Given enough divergence time, different deleterious mutations inherited from the ancestral population will reach different frequencies in each subpopulation. We investigate the genotype-phenotype relationship in the resulting mixed population. Each individual is characterized by the probability x that the two alleles it has at a given locus originate from the same subpopulation, depending on the number g of panmictic generations following the admixture. The quantity 1 - x can be interpreted as the individual "degree of outbreeding," i.e., the "disparity between the genome of the two parents" (COULSON et al. 1999 ).

Just after the admixture (g = 1) the correlation is due to the coexistence of two inbreeding classes, C_w and C_b, with equal probabilities 1/2. C_w (w stands for "within") is composed of individuals with both parents originating from the same subpopulation (x = 1) while C_b individuals (b stands for "between") have one parent in each of the two subpopulations (x = 0). C_w individuals are expected to have lower fitness than C_b individuals as they are more likely to be homozygous for their deleterious mutations. They also have more similar alleles at a neutral marker locus than C_b individuals.

Each generation after the contact, new inbreeding classes C_x appear as a consequence of recombination. At generation g after the contact, assuming free recombination for the sake of simplicity,

(11a)

(11b)

Let Q_w and Q_b be the probability of identity-in-state of genes from C_w and from C_b, respectively. f_x, the inbreeding coefficient of class x, is defined as

(12)

Neglecting mutation since admixture, Q_|x can be expressed as xQ_w + (1 - x) Q_b so that

(13)

where F_ST is the classical F-statistic (WRIGHT 1951 ; COCKERHAM and WEIR 1987 ). Note that F_ST is defined between the two subpopulations before the admixture and cannot be evaluated in samples of the current population (after the admixture). HFC studies are usually done within a single population and our model therefore concentrates on that situation. Although F_IS and F_ST are used in the context of the partial selfing and the admixture models, respectively, both models use the same definition of inbreeding, i.e., the increase in the probabilities of identity. By analogy with (3) the expected value of the fitness trait of an individual belonging to class x can be expressed as a function of f_x, or simply of x:

(14)

Here W₀ is the value of the fitness trait for C_b individuals (i.e., with x = 0) and the inbreeding load ß, the difference in fitness for C_b and C_w individuals, measures heterosis between the two subpopulations before the mixture. The random effect with mean 0 and variance ² is assumed independent of x. We now derive the correlation coefficient between heterozygosity and the fitness trait. Moments of the fitness trait W are derived on the basis of (11) and (14):

(15a)

(15b)

Moments of heterozygosity H are simply expressed as functions of probabilities of IIS:

(16a)

(16b)

The covariance term is calculated using the partition into inbreeding classes,

(17a) which yields

(17b)

The maximum correlation coefficient is

(18a)

At the first generation after the mixture (g = 1),

(18b)

The coefficient ²_max depends only on the ISS probabilities Q_b and Q_w (derived in Appendix B) but not on the inbreeding load, just as in the partial selfing model. In the presence of within-pedigree variance for the fitness trait, ² is as in (7b), replacing ²(f_k) by ²(x).

The model can be generalized to other cases, as for the partial selfing model. A summary of the analytical expressions for the squared correlation coefficients obtained under both population models is provided in Table 1.

Numerical parameters for the models:
Analytical results derived in the previous sections were explored numerically using Mathematica 3.0 programs (WOLFRAM 1996 ). The ranges of values explored for the various parameters of the model are presented below.

Mutation rates: To analyze the impact of mutation on heterozygosity-fitness correlations with respect to different mutation models, a large range of mutation rates is investigated, from 10^-9 to 10^-2 per locus and per generation. Mutation rates ranging from 10^-6 to 10^-2 are investigated to compare d² and H under the SMM, as this is the accepted range for a microsatellite locus. Assuming a substitution rate equal to 10^-9 per site per generation, total mutation rates ranging from 10^-9 to 10^-6 correspond to sequences of realistic length (few to several hundreds of base pairs). These mutation rates were thus considered when comparing p and H under the ISM.

Population parameters: Under the partial selfing model the standard situation modeled is a large population (n = 10³ individuals) with intermediate selfing rate (S = 0.4). The effect of a change in the population size (from 10³ to 10⁶) and the selfing rate (from 0 to 1) is analyzed. Under the admixture model the standard assumptions are small subpopulations (n = 10²) descended from a large ancestor population (N = 10⁴) after a long divergence time ( = 10⁴ generations). The effect of the subpopulation size (from 10² to 10⁴) and the joint influence of the ancestor population size and the divergence time in generations (assumed to be equal for the sake of simplicity and ranging from 10⁴ to 10⁶) are investigated. The inbreeding load does not need to be estimated here as we focus on maximum correlation coefficients, which do not depend on this parameter (see Table 1).

	RESULTS

TOP ABSTRACT THEORY RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C LITERATURE CITED

We first focus on the effects of the mutation models and mutation rates on the heterozygosity-fitness correlation for both population scenarios, evaluating the effect of population parameters. We then analyze the influence of genotypic indices on the genotype-fitness correlation obtained. Finally the form of the relation between the genotype and the fitness (the expected value of the fitness trait as a function of d²) is discussed under the SMM.

Impact of marker mutation on the heterozygosity-fitness correlation:
Under the partial selfing model, the maximal correlation coefficient increases with the marker mutation rate for all mutation models studied (Fig 1A). This is not surprising since this correlation is an increasing function of marker gene diversity [(1 - Q₁) in Equation 7a], which itself increases with the mutation rate irrespective of the mutation model (see Appendix A). For low mutation rates (<10^-5), the correlation coefficient is never higher than a few percent, in agreement with most experimental data (DAVID 1998 ). Mutation models have a low impact except for high mutation rates (>10^-5), which are realistic only for microsatellite markers. For high mutation rates the correlation coefficient increases with the value of K in the K-alleles model (KAM) and is maximal for the IAM. The SMM behaves like a KAM with large K (K > 10, data not shown).

View larger version (16K): In this window In a new window Download PPT slide   Figure 1. Influence of mutation on the correlation coefficient between heterozygosity at one marker locus and fitness. IAM, infinite alleles model of mutation; SMM, stepwise mutation model; 4AM, four-alleles model; 2AM, two-alleles model. The variance for the fitness trait () was assumed to be equal to 0 within inbreeding class. (A) Partial selfing model. The population size (n) is 103 and the selfing rate is 0.4. (B) Admixture model. Results are given at the first generation following the admixture. The subpopulation size (n) is 102, and the ancestor population size (N) and the divergence time () in generations are 104.

A different picture is obtained under the admixture model, as the maximum correlation coefficient exhibits a maximum for an intermediate mutation rate (Fig 1B). This "optimal mutation rate" is of the order of 1/ (the divergence time; Fig 1B and data not shown). The correlation coefficient obtained under the admixture model may be much higher than that obtained under the partial selfing model. Under the partial selfing model the correlation is limited by the fact that all genotypes, be they homozygotes or heterozygotes, are represented in all inbreeding classes although the frequency of heterozygotes is halved in each generation of selfing. In contrast, under the admixture model, parameter values can be found such that all homozygotes are in class C_w and all heterozygotes are in class C_b, which allows the maximum correlation (in the absence of environmental variance in fitness) to approach unity. However, assuming freely recombining loci, this maximum correlation is halved each generation following the admixture (Equation 18aEquation 18b). Furthermore, in natural populations, within-inbreeding class variance for the fitness trait may weaken the real correlation [Equation 7b, replacing ²(f_k) by ²(x)]. Mutation models rank in the same order as those under the partial selfing model with regard to the correlation coefficient, except with very high mutation rates (>10^-3) for which the SMM provides a stronger correlation than the IAM. Ultimately, for both population models studied, marker mutation influences HFC primarily through its effect on marker diversities: (1 - Q₁) for the partial selfing model (Equation 7a) and (1 - Q_w) and (1 - Q_b) for the admixture model (Equation 18aEquation 18b).

Impact of population parameters on the heterozygosity-fitness correlation:
Under the partial selfing model, the correlation coefficient increases with the population size and the selfing rate, until S is very close to 1 (Fig 2), whatever the mutation rate and model (data not shown). OHTA and COCKERHAM 1974 obtained qualitatively similar results for the effect of the selfing rate on associative overdominance in an infinite partially selfing population. However, CHARLESWORTH 1991 found that associative overdominance is maximal for intermediate selfing rates. This discrepancy probably results from the measure of associative overdominance chosen. CHARLESWORTH 1991 studied the apparent selection coefficient for heterozygotes at the neutral marker, which is expressed in phenotypic units and thus depends on the inbreeding load. The maximal heterozygosity-fitness correlation coefficient used here is independent of the inbreeding load ß, unlike the covariance term, which scales with ß (Equation 6b). Using the latter as a measure of associative overdominance, we obtain a maximum associative overdominance for intermediate selfing rates (data not shown), just as found by CHARLESWORTH 1991 . Taking into account within-pedigree variance for the fitness trait ( > 0 in Equation 4b) the correlation is weakened and becomes an increasing function of the inbreeding load (Equation 7b).

View larger version (15K): In this window In a new window Download PPT slide   Figure 2. Influence of population parameters on the correlation coefficient between heterozygosity at one marker locus and fitness, under the partial selfing model and IAM. ( = 0). The mutation rate is 10-6. The population size n varies from 103 to 106.

In the admixture model, the correlation coefficient decreases with the subpopulation size, particularly when the mutation rate is high (Table 2). The correlation increases when the ancestral population size and the divergence time increase simultaneously, particularly for low mutation rates. The effect of the inbreeding load ß is similar to that in the partial selfing scenario (Equation 7b).

Impact of the genotypic index on the correlation with the fitness trait:
The squared difference in repeat units vs. heterozygosity under the stepwise mutation model: Under the partial selfing model, for low mutation rates, d² and heterozygosity H are equally correlated with fitness (Fig 3A). Increasing the mutation rate strongly increases the correlation with heterozygosity but not with d² (Fig 3A), irrespective of population size or selfing rate (data not shown). This can be explained by the distribution of the expected value of the fitness trait conditioned on the value of d². Using (C3) and (C4) in Appendix C we found that all heterozygotes (nonnull d² values) correspond numerically to the same expected fitness (e.g., W = 0.89 with W₀ = 1, ß = 1, u = 10^-4, and other parameters as in Fig 3A), which is higher than that of homozygotes (d² = 0) (e.g., W = 0.72 with W₀ = 1, ß = 1, u = 10^-4, and other parameters as in Fig 3A). Therefore, heterozygosity is as informative as d², and it provides a stronger correlation due to its lower variance.

View larger version (16K): In this window In a new window Download PPT slide   Figure 3. Correlation coefficients between heterozygosity or d2 at one marker locus and fitness under the SMM ( = 0). (A) Partial selfing model. The parameters used are as in Fig 1A. (B) Admixture model. The parameters used are as in Fig 1B except that n = 103.

Under the admixture model, however, the correlation between heterozygosity and fitness is weakened for high mutation rates (Fig 1B and Fig 3B) whereas the correlation with d² is almost insensitive to the mutation rate (Fig 3B). With a small subpopulation size (n = 100) and mutation rate (u = 10^-4), the correlation between the fitness trait and d² is weaker than that with heterozygosity, as in the partial selfing model (data not shown). In contrast, when the subpopulation size and the mutation rate are large enough (so that nu >> 1) d² provides a better correlation than H (Fig 3B). Again, this can be explained by the distribution of the expected value of the fitness trait conditional on d², derived according to Appendix C (Fig 4). Indeed in the first case (nu << 1), d² has little value for predicting fitness once it is >1. However, when nu >> 1, the expected value of fitness increases progressively with the value of d², so that in this case d² is far more informative than H.

View larger version (10K): In this window In a new window Download PPT slide   Figure 4. The expected value of the fitness trait as a function of d2 value under the admixture model assuming the SMM ( = 0). The parameters used are as in Fig 3B.

The number of nucleotide differences vs. heterozygosity under the infinite-sites model: Under the partial selfing model, given the low mutation rate assumed for a neutral sequence, the same correlation with fitness is obtained whether heterozygosity H or the number of nucleotide differences p is considered (Fig 5A). Only very high mutation rates (>10^-4), corresponding to unreasonably long sequences, would introduce a difference, with a higher correlation with H than with p (data not shown).

View larger version (15K): In this window In a new window Download PPT slide   Figure 5. Correlation coefficients between heterozygosity or the number of substitutions and fitness under the ISM ( = 0). (A) Partial selfing model. The parameters used are as in Fig 1A. (B) Admixture model. The parameters used are as in Fig 1B except that n = 104.

Under the population mixture model, identical correlation coefficients are obtained with p and H except with large subpopulation size (>10³) and large mutation rates (>10^-6; data not shown). In this case p is slightly less correlated to fitness than H (Fig 5B). Again, only very high mutation rates (>10^-4) would result in a higher correlation with p than with H (data not shown). Unsurprisingly, these results are qualitatively similar to that obtained in the comparison of d² and H under the SMM with equivalent mutation rates.

The expected value of the fitness trait as a function of d²:
The expected shape of the genotype-phenotype relationship is not linear (see Fig 4), with fitness plateauing for high d²-values. This indicates that linear regression models previously used in empirical studies are not appropriate. We therefore suggest using a nonlinear regression based on the logistic equation

(19)

where j is the square root of the value taken by d², W() is the maximum value of W, asymptotically reached for very high d² values, W(0) is the value associated with d² = 0, and r measures the rate at which the fitness trait value saturates as a function of d². This can be rewritten:

(20)

For the partial selfing model, (d²) can be interpreted as a microsatellite-based indirect measure of inbreeding depression. For the admixture model it is interpreted as an indirect measure of heterosis between the two subpopulations before the mixture. Numerical data obtained using Appendix C were fitted to this nonlinear model, for a large range of parameters defined in Table 3 and Table 4 (data not shown). This model explains 95–100% of the variance in fitness obtained for both population models studied (neglecting within-inbreeding class variance for the fitness trait). It thus provides a satisfactory description of the genotype-phenotype relationship on the basis of three parameters. Under the partial selfing model, (d²) increases with the marker mutation rate, the selfing rate, and the population size (Table 3). Under the mixture model, (d²) increases with the marker mutation rate and when the divergence time and the population size increase simultaneously. It decreases with the subpopulation size (Table 4).

View this table: In this window In a new window   Table 3. (d2) under the partial selfing model

View this table: In this window In a new window   Table 4. (d2) under the admixture model

	DISCUSSION

TOP ABSTRACT THEORY RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C LITERATURE CITED

A unified framework for studying genotype-fitness relationships:
To our knowledge, mutational processes at marker loci have not previously been incorporated into analytical models of associative overdominance, although they are recognized as relevant for empirical issues (PEMBERTON et al. 1999 ). Evaluating heterozygosity-fitness correlations requires the evaluation of probabilities of identity-in-state, which depend both on the demographic scenario assumed and on mutational processes and have been extensively used to address population structure issues (reviewed in ROUSSET 2001 ). This could be done in a unified framework using the same formalism as previously used (ROUSSET 1996 ). It is equivalent to a formalism based on the computation of distributions of coalescence times but does not require an explicit computation of such distributions. Under this formalism, neglecting within-inbreeding class variance in fitness for a given fitness model, the maximal correlation coefficient between a fitness trait and a neutral genotype does not depend on the inbreeding load (the reduction in fitness associated with complete inbreeding), as already reported for other demographic scenarios (BIERNE et al. 2000B ). There is therefore no need for an explicit expression of the inbreeding load, which depends on the genetic architecture of the fitness trait and on the population scenario. Rather than exploring the continuum of possible historical scenarios exhaustively, we chose to exemplify a short-term (with only very recent coalescence events) and a long-term (with coalescence events deeper in the pedigree of individuals) inbreeding process with, respectively, the partial selfing and the admixture model.

The fitness model we used (Equation 3 and Equation 14) was first proposed by MORTON et al. 1956 . Neglecting purging selection and disequilibrium between selected loci and assuming additive effects among fitness loci (or multiplicative effects, considering log-fitness rather than fitness), we obtained the Morton model for fitness under various genetic architectures and various population scenarios (CHARLESWORTH and CHARLESWORTH 1987 , CHARLESWORTH and CHARLESWORTH 1999 ; BIERNE et al. 2000B ; WHITLOCK et al. 2000 ).

Variance for the fitness trait may exist within an inbreeding class, particularly when recombination is limited, due to the segregation of chunks of chromosomes rather than independent loci. In natural populations, environmental effects may also increase this variance. If this is large enough compared to genetical effects, HFC would be weakened and also would depend on the inbreeding load, which should then be estimated. Empirically estimating within- and between-inbreeding class variance for the fitness trait might be difficult in natural populations, but would be particularly useful to assess the validity of our model.

Choosing appropriate marker genes:
Microsatellites have recently been used as a tool to infer fitness differences due to variation in the level of inbreeding between individuals (PEMBERTON et al. 1999 ). Are microsatellites good markers to address such issues?

We have shown that HFC due to partial selfing increases with the marker mutation rate (and thus with marker variability), making microsatellites good markers in such a situation. Partial selfing is an example of short-term inbreeding generated by the current mating system. One generation of outcrossing reduces individual inbreeding to zero, so that only selfing since the most recent outcrossing event (i.e., a few generations) needs to be taken into account. However, short-term inbreeding can also be due to recent demographic episodes such as bottlenecks. A sudden reduction in population size produces random inbreeding among individuals, i.e., mating between kin due to chance, as opposed to systematic inbreeding due to the mating system (MALECOT 1969 ). BIERNE et al. 2000B analyzed theoretically the situation of a population experiencing a recent and drastic bottleneck sustained over a few generations. They found strong HFC, but it is transient due to the rapid homogenization of the inbreeding level among individuals after a few generations of a sustained bottleneck. BIERNE et al. 2000B did not use an explicit mutation model for the marker genes. However, they found that highly variable markers produce larger HFC than less variable ones, consistent with our result under partial selfing. Heterozygosity-fitness correlations increase with marker diversity and thus with marker mutation rates, whenever they are related to short-term inbreeding. Microsatellite loci seem therefore ideal to investigate fitness consequences of short-term inbreeding.

However, when the origin of HFC is long-term inbreeding (e.g., admixture of anciently diverged small subpopulations), the correlation increases up to a certain mutation rate but then decreases (Fig 1B). Using highly mutable markers can thus lead to confusing results under a long-term inbreeding scenario. The explanation is that even highly "inbred" genotypes can have heterozygous marker loci if marker mutation rates are sufficiently high. Marker heterozygosity is therefore not a good index of individual fitness. There is thus an "optimal mutation rate," which decreases as the timescale of the population inbreeding history increases. In conclusion, one should use molecular markers whose mutational process scales roughly with the assumed cause of inbreeding. However, empirical studies provide only imprecise estimates of mutation rates and a temporal scale of inbreeding. In practice we suggest only that empirical studies should avoid using markers that mutate at an obviously inappropriate rate.

Mutation models of marker genes play little role in HFC, except when the mutation rate is high. In the infinite-alleles model there is no homoplasy, whereas the K-alleles model is characterized by increasing homoplasy with decreasing K-values. The comparison between IAM and KAM, and among various K-values, suggests that homoplasy is a limiting factor for HFC, irrespective of the population scenario. Homoplasy represents a loss of information about identity-by-descent, making marker genotypes less closely related to the inbreeding level of individuals. As the homoplasy generated by the stepwise mutation model is different from that produced by a KAM, the results obtained for these two mutation models cannot be directly compared.

Choosing an appropriate genotypic index for microsatellite data:
It has been suggested recently that d² may be a more appropriate index than heterozygosity when analyzing fitness consequences of inbreeding (COULSON et al. 1999 ; PEMBERTON et al. 1999 ). Heterozygosity should be adequate to distinguish inbred and outbred individuals, whereas d² may give the additional opportunity to detect "highly outbred" individuals having high d² values, i.e., large coalescence time for their two microsatellite alleles, vs. "moderately outbred" individuals (COULSON et al. 1999 ).

In a red deer population from the Isle of Rum, COULSON et al. 1998 found a positive correlation between birth weight and mean d², although not between birth weight and heterozygosity. They argued that this could be due to the demographic history of the population, namely population divergence followed by mixing (our admixture model), which is supported by historical data. However, we have shown that, in our admixture model, at least, the conditions under which d² would lead to a better correlation than H in the admixture model are restrictive (a very high marker mutation rate and large subpopulation sizes and divergence time) and furthermore they do not seem to apply to what is known of the red deer population of the Isle of Rum (PEMBERTON et al. 1999 ). Divergence times of 100 generations, with large subpopulation size, parameters that appear realistic for such a population, always lead to better predictive power of H as compared to d² (not shown). Note that the fact that there are only two classes of inbreeding and fitness in our admixture scenario just after the contact does not necessarily imply that a binary variable such as H will perform better than a quantitative variable such as d². If the proportion of high fitness individuals (class C_b) within a class of d² increases progressively with d² (Fig 4, with u = 10^-2) rather than in a stepwise way (in most cases), d² performs better than H.

One can suggest another explanation for the stronger correlation with fitness obtained using d² as compared to H. First, we cannot exclude the possibility that a more sophisticated demographic scenario may generate this result under less restrictive conditions than the admixture scenario formulated here, although we do not have a precise idea of what such a scenario should be. Alternatively, the results of COULSON et al. 1998 , COULSON et al. 1999 may be due to low homozygosity in their population, making heterozygosity an inadequate genotypic index when compared to d². Finally these results could simply be due to chance, i.e., random variation. Recently, several studies have found fitness-related traits, such as adult breeding success in the red deer (SLATE et al. 2000 ) or survival and parasite resistance in the Soay sheep (COLTMAN et al. 1999 ), to be more strongly correlated with heterozygosity than with d². In a captive wolf population with known inbreeding levels, HEDRICK et al. 2001 found that d² was less predictive of the known inbreeding coefficient than microsatellite heterozygosity. This supports the hypothesis that d² is generally not more powerful than heterozygosity to detect fitness consequences of inbreeding.

It seems generally unlikely that d² will provide a higher correlation with fitness than H. For low mutation rates, heterozygosity and d² are equivalent under both our population models, and for high mutation rates, heterozygosity always provides a higher correlation than d² with fitness under the partial selfing model. Under the admixture model only, d² is more correlated with a fitness trait than heterozygosity under restrictive conditions: large mutation rate, divergence time, and subpopulation size. These conditions are expected to be associated with a relatively low inbreeding load (ß) in the mixed population, because there is a small probability of fixation, or of sufficient change in frequency, of a deleterious mutation in a large subpopulation. The inbreeding load does not influence the maximum correlation (Equation 18a) but affects the actual correlation, taking into account within-inbreeding class variance for the fitness trait. Given that sources of variation in fitness other than inbreeding must exist in natural populations (generating in Equation 4b and Equation 15b), this small inbreeding load weakens the actual correlation (see Equation 7b), thus reducing the chance of detecting it.

Finally, a conclusion of our analysis is that there seems to be little theoretical reason to use d² instead of H when analyzing correlations between microsatellite genotype and fitness, even for long-term inbreeding scenarios. When doing so, however, one should choose an appropriate regression model to analyze genotype-fitness correlations (Equation 19 and Equation 20) rather than the linear model used in previous studies.

	ACKNOWLEDGMENTS

We thank Nicolas Bierne, Sylvain Glémin, Philippe Jarne, Sally Otto, Josephine Pemberton, John Slate, and John Thompson for stimulating discussions and helpful comments on a previous draft of the manuscript. This work was supported by funds from the Centre National de la Recherche Scientifique (CNRS; Jeune Equipe program) to P. Jarne.

Manuscript received May 23, 2001; Accepted for publication October 1, 2001.

	APPENDIX A

TOP ABSTRACT THEORY RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C LITERATURE CITED

IDENTITY-IN-STATE IN THE PARTIAL SELFING MODEL
We derive probabilities of identity-in-state within and between individuals under the infinite sites, the stepwise mutation, and the K-alleles mutational models (see ROUSSET 1996 ). Generating functions of allele size difference (under the SMM) or of the number of nucleotide differences (under the ISM) are also derived and used to obtain the moments of d² and of p.

IAM and ISM:
Recursions for the IIS probabilities under the IAM or the ISM are

(A1)

with = (1 - u)². Note that in the ISM or the IAM, IIS and identity-by-descent (IBD) are equivalent. Q₀ and Q₁ are the equilibrium values of the system above:

(A2)

Under the ISM the generating function ₀ (resp. ₁) of the probability p_j,0 (resp. p_j,1) that two randomly chosen alleles within an individual (resp. between individuals) differ by j nucleotide sites is defined as ₀(z) = ⁺_j=0p_j,0z^j (resp. ₁(z) = ⁺_j=0p_j,1z^j). The effect of mutation is to change both generating functions (z) by a factor: r(z) = (1 - u + uz)². We also introduce _0|k, which is the value of ₀ conditional on inbreeding class k.

The recursions on are then mathematically similar to those for IBD (ROUSSET 1996 ):

(A3)

₀ and ₁ are the equilibrium points of the system above, and neglecting mutation since the last outcrossing event _0|k is easily derived:

(A4)

By differentiating these generating functions one obtains the various moments of p:

(A5)

SMM:
To derive Q₀ and Q₁ under the SMM it is useful to introduce ₀ (resp. ₁), the generating function of the probability p_j,0 (resp. p_j,1), that two randomly chosen alleles within an individual (resp. between individuals) differ by j steps. Because steps can take either nonnegative or negative values, ₀(z) = ⁺_j=-p_j,0z^j (resp. ₁(z) = ⁺_j=-p_j,1z^j). We also define _0|k, the value of ₀ conditional on class k. The effect of mutation is to change the generating functions by a factor r(z) = (1 - u + ())², and ₀ and ₁ are calculated as in (A4) using this new r(z). By definition, we have p_0,0 = Q₀ and p_0,1 = Q₁. Furthermore, p_-j,0 = p_+j,0 = Pr(d² = j²) for all j 0. Therefore,

and neglecting mutation since the last outcrossing event,

We use the inverse Fourier transform in j of a given function f:

Q₀ and Q₁ can computed as

(A6)

L_j(f) are numerically calculated with Mathematica 3.0 (WOLFRAM 1996 ).

By differentiating the generating functions one obtains the various moments of d²:

(A7)

KAM:
The recursions for IIS probabilities are

(A8)

in which = (1 - u) + , Q^'_j,t = 'Q_j,t + , and ' = ( )².

The equilibrium values of the system are

(A9)

	APPENDIX B

TOP ABSTRACT THEORY RESULTS DISCUSSION APPENDIX A APPENDIX B APPENDIX C LITERATURE CITED

IDENTITY-IN-STATE IN THE MIXTURE MODEL
We derive probabilities of identity-in-state under the ISM, SMM, and KAM (see ROUSSET 1996 ) focusing on C_w, C_b, and C_x individuals. Generating functions of allele size difference (under the SMM) or of the number of nucleotide differences (under the ISM) are also derived, as they are useful to derive moments of d² and of p. In contrast to the partial selfing model there is no need to distinguish identity within (subscript 0) and between (subscript 1) individuals, as random mating is assumed.

IAM and ISM:
The probability of identity (IIS and IBD are confounded here) for two alleles randomly drawn from the ancestral random mating population of size N at inbreeding equilibrium is Q_a with = (1 - u)². Let two subpopulations from the ancestral population diverge without gene flow for generations. Just before the admixture, the probabilities of identity of two randomly chosen alleles are Q_w() = Q_eq + (())(Q_a - Q_eq) within a subpopulation and Q_b() = Q_a between subpopulations, with Q_eq = . After the admixture, neglecting mutations since the admixture event and considering an infinite mixed population for the sake of simplicity,

(B1)

Under the ISM we obtain analogous expressions for generating functions, as we did in APPENDIX A. Let _w (resp. _b) be the generating function of the probability p_j,w (resp. p_j,b) that two randomly chosen alleles that originated from the same subpopulation (respectively from the two subpopulations) differ by j nucleotide substitutions [_w(z) = ⁺_j=0p_j,wz^j]. Recursions on (z) are mathematically identical to those for IBD, taking into account mutation through the factor r(z) = (1 - u + uz)². Therefore, just before the admixture we obtain

within subpopulations and _b,(z) = (r(z))_a(z) between subpopulations, with