Introduction

Self-fertilization is the fusion of gametes from a single genetic individual. It is common in angiosperms, a group with <10% dioecious species (Renner and Ricklefs, 1995), as well as in other plants, such as mosses (Eppley et al., 2007), and in hermaphroditic animals, such as mollusks or trematodes (Jarne and Auld, 2006). It also occurs in fungi. Estimating selfing rate is an important issue in population and evolutionary biology for two reasons. First, selfing substantially affects the distribution of genetic variation within and among populations, and therefore the response of populations to selection (review in Jarne, 1995; Charlesworth, 2003). Second, the evolution of the selfing rate itself is a topic of central importance (review in Jarne and Charlesworth, 1993; Uyenoyama et al., 1993; Barrett, 2003; Goodwillie et al., 2005). Precise estimates are required to understand the evolution of selfing in relation to reproductive traits (for example, floral traits, copulatory behavior) and inbreeding depression, as well as to test a central prediction of genetic models of mating system evolution, that distributions of selfing rates across species should be U-shaped (Jarne and Auld, 2006, in animals; Goodwillie et al., 2005, in plants).

The interest in self-fertilization, especially in plants, goes back as far as the mid-nineteenth century (Darwin 1876). However, quantifying the selfing rate became possible only much later. The first attempt might be credited to Jones (1916) who compared parents and offspring genotypes at a dominant morphological marker in tomato, and provided a maximum estimate of the selfing rate. This might be considered as the ancestor of one of the two approaches used today, the progeny-array approach (PAA). The other main approach is based on the analysis of population genetic structure (population structure approach; PSA in what follows). The most well-known PSA method relies on the correspondence between the inbreeding coefficient and the selfing rate (Li, 1955). However, this method became popular only when estimating the inbreeding coefficient was possible, that is with the rise of molecular markers. Curiously, although the selfing rate has been estimated in an extremely large number of studies, these approaches have apparently not been reviewed since the late 1990s (Brown et al., 1989; Brown, 1990).

Our goal here is to fill this gap, that is to review the methods currently used for estimating the selfing rate in populations based on molecular markers, covering progresses in both molecular and statistical methods. After briefly outlining the available markers, we present how they can be used in PSA and PAA. We mention their strengths and weaknesses, including methodological and technical issues, and argue that they can fruitfully be run together. Note that estimating the selfing rate is often embedded in a wider context, for example in surveys of among-population genetic structure or paternity analysis (see for example, Blouin, 2003), but we restrict our comments to the two approaches mentioned above.

Molecular markers

Mating system parameters, such as the selfing rate, have been estimated using morphological markers with simple genetic basis, especially before the rise of molecular markers. They include flower color in plants (Fryxell, 1957; Cruzan, 1998) or mantle pigmentation in freshwater snails (Vianey-Liaud, 1997). However, selection acts on these markers, often in a rather complex way (for example, Irwin and Strauss, 2005), and it is difficult to get more than one marker per species. The most popular tools are now molecular markers that are assumed to be selectively neutral (see Avise, 2000). The choice of an appropriate marker depends on the same broad criteria whatever the method used (Supplementary Appendix 1): (1) Co-dominant markers are preferable, because it is possible to unambiguously identify heterozygous individuals. Dominant markers can however be used with the PAA (but not with the PSA). (2) Mendelian transmission across generations is required, and is usually assumed for molecular markers (for example, Selkoe and Toonen, 2006). Deviations from co-dominance and Mendelian segregation can in principle be evaluated using progeny-arrays. (3) The precision of all estimates increases with genetic diversity, meaning that highly polymorphic markers are preferred. This requirement is often difficult to meet in highly selfing species, which are usually less polymorphic than their outbred counterparts (Jarne, 1995; Charlesworth, 2003).

The three most commonly used markers are allozymes, microsatellites and AFLP (amplified fragment-length polymorphism), and relevant characteristics are highlighted in Supplementary Appendix 1. Allozymes remain by far the most widely used, totaling 82% of entries in the database used by Goodwillie et al. (2005), and 97% of that in Jarne and Auld (2006). Based on literature searches (using ISI Web of Knowledge) and Goodwillie et al. (2005), about 15–25 estimates of selfing rate in plants have been published yearly since 1990 using allozymes. Allozymes are co-dominant and five to ten polymorphic loci can easily be scored. However, their generally limited polymorphism in highly inbreeding species is a major incentive for developing microsatellites (for example, Viard et al., 1996). Microsatellites exhibit the same general characteristics of allozymes, but are usually more polymorphic (Jarne and Lagoda, 1996; Estoup and Angers, 1998). Although their use remains limited, with about 6% of the 469 entries in Goodwillie et al. (2005), it is increasing with about 10 papers per year since 2001. Although 10 loci can easily be characterized through molecular cloning, this is technically more demanding and time consuming than developing allozymes. AFLPs (Vos et al., 1995), on the other hand, do not require cloning. They are considered as dominant markers with two alleles per locus (present/absent), although the probability that different bands actually represent alleles of the same locus is generally unknown. These limitations might, to a certain extent, be counterbalanced by the large number of loci scored (several tens to hundreds). Dominance imposes PAA as the only possible method for estimating the selfing rate with AFLPs. Very few studies have been published to date (for example, Thompson and Ritland, 2006).

A major problem with molecular markers is technical bias (see Supplementary Appendix 1). Selfing rate analysis assumes strict correspondence between a scored phenotype (for example, a band on a gel) and an allele. Unfortunately even cross-checked data sets are not error free (Hoffman and Amos, 2005; Pompanon et al., 2005). Most technical biases can be seen as problems of partial dominance, whereby heterozygotes are read as homozygotes for one of the two alleles. This is the case for null alleles that affect both allozymes (for example, David et al., 1997) and microsatellites (see Dakin and Avise, 2004; Selkoe and Toonen, 2006). The same effect is produced whenever two bands in a heterozygote are too close to be visually separated, which can occur in the case of band stuttering at microsatellite loci, or when microsatellite alleles compete for amplification (short allele dominance; Selkoe and Toonen, 2006). In addition, even if heterozygosity is correctly detected, erroneous allele identification affects the expected genotypic frequencies, hence the estimate of selfing rate, whatever the approach. For example if allele A is occasionally read as a new allele B, the observed lack of AB heterozygotes will automatically increase the apparent selfing rate. Other technical problems are related to pattern repeatability and Mendelian transmission across generations (for example, alleles scored in offspring but not in parents). Such problems are probably more acute with AFLPs, but should not be underestimated in the case of microsatellites and especially allozymes (for example, in relation to post-translational modifications of proteins). The influence of technical biases on estimates of the selfing rate is considered in more detail below.

On the whole, microsatellites, arguably as technically robust as allozymes, but usually more polymorphic, are currently the most appropriate markers for estimating inbreeding, and their main limitation is their technical and economic cost.

Population structure approach

The population structure approach is based on samples taken directly in natural populations; all individuals sampled are taken from the same generation. We first describe the classical analysis of such data using single-locus inbreeding coefficients, and then more recent multilocus techniques.

Using single-locus inbreeding coefficients

A simple relationship connects selfing rate, S, to the inbreeding coefficient, F (Fyfe and Bailey, 1951; Li, 1955; Charlesworth, 2003):

This equation assumes inbreeding equilibrium in an infinite population of diploid organisms in which selfing is the sole source of inbreeding. Outcrossing occurs at random, and there is no selection, especially no inbreeding depression. A similar formula is available for tetraploids (Ronfort et al., 1998) and situations in which both gametophytic and sporophytic selfing may occur, such as in mosses (see Eppley et al., 2007). A strong advantage of this approach is that F is almost systematically estimated in analyses of population genetic structure, based on co-dominant markers, providing an ample source of data. For example, estimates of selfing in animals are essentially derived from this approach (Jarne and Auld, 2006).

Various methods and softwares are available for estimating F using genotypes at single loci (Ayres and Balding, 1998; Excoffier, 2001; Rousset, 2001), the most popular estimate being that of Weir and Cockerham (1984). Resampling techniques can be used to assess whether F, and hence the selfing rate, differs from 0 (Weir and Cockerham, 1984; Goudet, 1995). The single-locus sampling variance of F̂ is a function of F, n (the number of individuals sampled) and allelic frequencies (Curie-Cohen, 1982; see Supplementary Appendix 2). The variance of Ŝ follows:

Var(Ŝ) is therefore of the order of 4Var(F̂) when F is close to 0, and close to Var(F̂)/4 when F=1. The essential points to remember are that the error on Ŝ decreases with n and He (gene diversity), and generally with S itself. It can be shown that increasing n is more effective than increasing the number of loci for improving the precision in Ŝ for the same genotyping effort (Supplementary Appendix 2). However, using several unlinked loci is always a good idea because it limits the impact of locus-specific biases (that is nonneutrality or technical problems). Comparing estimates based on different loci provides a simple evaluation of the overall robustness of the method. Variances and confidence intervals across loci can also be obtained (Excoffier, 2001). For example, Weir and Cockerham (1984) proposed using the jackknife. Finally, when possible, loci with large values of He should be typed in priority (Supplementary Appendix 2).

Multilocus approaches

The inbreeding coefficient method can combine information from several loci as a weighted average. This approach certainly increases the precision of F estimates, but genotypic associations among loci are not exploited. Such multilocus information is, however, used by more recent methods (see David et al., 2007). In principle, the simplest way to incorporate this information is to find the selfing rate maximizing the likelihood of all observed multilocus genotypes. In practice however, many variables including allelic frequencies at all loci must be estimated as well. A useful shortcut is to evaluate allelic frequencies first, and then to find S assuming that they are known exactly. This is the approach followed, for example by Enjalbert and David (2000), for estimating S in artificial wheat populations whose founders have been genotyped. Criscione and Blouin (2006) proposed a different approach based on individual inbreeding coefficients estimated from actual population samples. The expected distribution of these coefficients in outcrossed and selfed offspring can be simulated and used to estimate the maximum likelihood function of S conditional on the observed distribution of inbreeding coefficients. However, the method allows at most a single generation of selfing in individual pedigrees, and is therefore not appropriate for high selfing rates.

David et al. (2007) proposed a novel method based on the distribution of multilocus heterozygosity that does not require joint estimation of allelic frequencies. Inbreeding generates not only heterozygote deficiencies, but also identity disequilibria, that is correlations of heterozygosity across pairs of loci (Weir and Cockerham, 1973). In other words, inbred individuals tend to be less heterozygous at all loci than outbred ones. If we consider two loci, the frequency of doubly heterozygous genotypes (and that of doubly homozygous genotypes) is higher than expected under independent assortment (Weir and Cockerham, 1973; David et al., 2007):

where h1 and h2 are indicators of heterozygosity at the two loci, E() stands for the expectation and g2, the second-order heterozygosity disequilibrium (or identity disequilibrium), gives the relative excess of genotypes heterozygous at two loci. David et al. (2007) generalized this idea to several loci and derived a multilocus estimator of S assuming inbreeding and linkage equilibrium. The equation

can be used to derive S from g2, and g2 can be estimated from the distribution of multilocus heterozygosity at a set of loci (David et al., 2007). Note that the only extra conditions required to derive (4), compared to (1), is linkage equilibrium. The sampling variance of Ŝ is comparable to that obtained using inbreeding coefficients and decreases when the number and genetic diversity of loci increase. David et al. (2007) also provide a maximum-likelihood estimate of S, based on the same general idea. This method is implemented in the software RMES (available at ftp://ftp.cefe.cnrs.fr/).

Estimates of the selfing rate can also be derived from the long-term effect of selfing on effective recombination and hence linkage disequilibria. For example, Cutter (2006) proposed an approach (detailed in Supplementary Appendix 5) using sequence data, assuming both drift/recombination and inbreeding equilibrium, and requiring some knowledge of recombination rates and effective population size. Given these assumptions, this method provides a long-term estimate of S. The number of assumptions and input parameters restrict the use of such methods and make them unlikely to produce precise estimates, especially when population size is limited, except perhaps for low selfing rates (see Supplementary Appendix 4).

Progeny-array approach

This approach is based on the comparison of mother and offspring genotypes. Ignoring mutation, A1A2 offspring with an A1A1 mother must be outcrossed and their frequency can provide a minimal estimate of the outcrossing rate (1−S), which of course depends on allelic frequencies. This is the essence of the PAA model for a single locus, set out in Supplementary Appendix 6. This model, or variants thereof, has been used to study plant mating systems since the premolecular era (Fryxell, 1957). Formal maximum likelihood estimation was initiated by Fyfe and Bailey (1951), and later extended to co-dominant loci (allozymes) by RW Allard and associates (for example, Brown and Allard, 1970; Clegg et al., 1978). A key insight was to incorporate multilocus information (Ritland and Jain, 1981; Shaw et al., 1981), and a good part of subsequent developments of the maximum likelihood method is due to the continuous efforts of Ritland (Ritland, 1986, 1989, 2002; Thompson and Ritland, 2006), and others (for example, Schoen and Clegg, 1986). This was paralleled by a marked increase in the number of studies estimating the selfing rate in natural plant populations (see Schemske and Lande, 1985; Goodwillie et al., 2005). However, the PAA began to be used a few years ago only in animals (see Jarne and Auld, 2006).

The PAA relies on several assumptions, including absence of population substructure and of selection between fertilization and sampling (Table 1; Supplementary Appendix 6). The main difficulty in estimating S is that a nonnegligible proportion (say, β) of outcrossed individuals have genotypes compatible with selfing (Table 2). This proportion can be inferred (using estimated allele frequencies) and used to correct estimates of outcrossing rates, as for example proposed by Shaw et al. (1981):

with n nonambiguously outcrossed offspring out of N in the array. This simple equation illustrates the advantage of using multilocus information. Large numbers of loci increase the probability of detecting outcrossing, hence increasing n and decreasing β, reduce the uncertainty on S. Such multilocus estimators rely on the same assumptions than the single-locus ones, plus the absence of linkage disequilibria. Following the same logic, Ritland and Jain (1981) proposed a more general, likelihood-based framework for the PAA, which has subsequently been refined (Ritland, 1986, 1989, 2002) (see Supplementary Appendix 6). The model was embedded in a more general framework and some assumptions were relaxed, allowing for a wider set of biological situations and introducing new parameters:

  • The initial model was intended for diploids assessed with co-dominant markers in a mixed-mating model. Extensions were developed for tetraploids with tetrasomic inheritance and no double reduction at meiosis (Murawski et al., 1994; Thompson and Ritland, 2006), automixis and apomixis (Thompson and Ritland, 2006), and dominant markers (Ritland, 1990b; Thompson and Ritland, 2006).

  • S can be estimated at different sampling levels, for example individuals and/or families (Ritland and El Kassaby, 1985). This opens the possibility of studying stratified populations, for example the case when two sexual morphs might exhibit different selfing rates, such as in gynodioecious plants or aphallic snails (Ritland, 2002).

  • It was already clear in Shaw et al. (1981) that differences between single-locus and multilocus estimates provide an estimate of the amount of biparental inbreeding. However, a more formal analysis was proposed through the effective selfing model (Ritland, 1984, 1986, 2002), an effective selfing event occurring either through uniparental selfing, or biparental inbreeding. In more genetical terms, effective selfing accounts for the co-ancestry between male and female contributions in outcrossed offspring. This was further elaborated in connection with the correlated-matings model (see below).

  • In the basic model, the mother genotype is compared separately to that of each offspring. Further information can be retrieved from the comparison of offspring genotypes, and this is the cornerstone of the correlated-matings model (Ritland, 1989, 2002). The selfing rate is here estimated together with three other parameters, namely rs (correlation of selfing), rp (correlation of outcrossed paternity) and paternal F. The correlation of selfing between two sibs is the covariance of selfing divided by S(1–S). rs can be interpreted as the fraction of sib pairs that are fully selfed, and it tends toward 1 when sibs are all selfed, or all outcrossed within a family. The correlation of outcrossed paternity is the proportion of full sibs among outcrossed sibs. The correlation rp is simply related to F and the probability of identity by descent for the two paternally derived gametes in the two progenies (Ritland, 1989). The inverse of rp can be interpreted as the effective number of fathers. A multilocus version of the correlated-matings model has been proposed (Ritland, 2002). One motivation is that the difference between single and multilocus estimates can be used to calculate biparental inbreeding, yet the uncorrected parameter will be an underestimate by a quantity which is a function of the number of loci (among other parameters), potentially biasing comparisons across studies.

Table 1 Influence of technical problems, biological processes (population structure and history, selection), marker characteristics and sample size on estimates of the selfing rate (and when relevant, on their variance) when using the single-generation single-locus (F), single-generation multilocus or progeny-array approaches
Full size table
Table 2 Genotypes of a mother and three progenies at four diploid loci (A–D) with up to three alleles
Full size table

One limitation of the PAA is that it is not always possible to obtain progeny-arrays from natural environments (see Supplementary Appendix 6). This is one of the reasons why it has rarely been used in animals (10 entries only in Jarne and Auld, 2006). On the other hand, its popularity among botanists (about 80% of entries in the database of Goodwillie et al., 2005) might be explained by several outstanding features, some of which are not a shared with the PSA. (1) The PAA provides estimates of S that are closer to primary selfing rates because offspring are usually sampled at juvenile stages. (2) Estimates of S can be stratified at lower levels, including families and individuals. However, when the variance in S among families is large, a large number of families should be studied in order to get a precise mean estimate of S. This seems to be about the case: we haphazardly extracted 30 studies, published from 1993 on, from Goodwillie et al. (2005). Out of the 32 species studied, the average number of families per population was 25.9 (s.d. 16.91), not much less than the ca. 30 individuals generally studied in PSA. In the 30 PAA studies, 17.1 (8.92) offspring were analyzed on average per family (the number was not correlated with the number of families in the study). (3) The PAA provides instantaneous estimates of S, which can be used to test for selective pressures acting at some point in time, whereas PSA estimates integrate various sources of inbreeding over several generations. (4) The comparison between single and multilocus estimates of S in the PAA allows detection of nonselfing sources of inbreeding (Ritland and Jain, 1981; Ritland, 1986, 2002). Note though that large sample sizes are required to obtain reasonable precision. A further reason of the success of the PAA is that the selfing rate, and correlated-matings parameters, can be estimated using the software MLTR (Ritland, 2002; http://www.genetics.forestry.ubc.ca/ritland/). The most recent version has been improved to accept large numbers of alleles per locus, a desirable quality when using microsatellites.

Technical and methodological biases

The main methods presented above (the inbreeding-coefficient approach, the multilocus approach of David et al. (2007) and the PAA) are not free of bias, and we distinguish technical biases affecting the chosen markers from methodological biases due to violation of model assumptions.

Technical bias

The technical problems affecting molecular markers (see above) produce biased estimates of the selfing rate, especially using inbreeding coefficients and PAA, as these approaches assume that genotypes are known without errors. Null alleles, short allele dominance and band stuttering generate heterozygote deficiencies, increasing F and apparent selfing (Table 1). Misreading alleles has the same effect, not because the observed heterozygosity decreases, but because the expected heterozygosity usually increases. On the other hand, the genotyping errors mentioned above should on average lead to underestimating S in the PAA. The reason is that a random error is more likely to produce a mother–offspring comparison compatible with outcrossing than with selfing.

The influence of technical bias on selfing rates derived from inbreeding coefficients can be quantified in the case of null alleles. The occurrence of such alleles is a more serious problem in outbred populations, because they remain undetected as null homozygotes are rare. For example, if the frequency of a null allele is 0.1 in an outbred population (S=0), the expected frequency of null homozygotes is 0.01, and they go unnoticed given the usual sample size of about 30–40 individuals per population. From the simple model presented in Supplementary Appendix 4, it can be shown that the expected difference between the estimated F̂ and the actual value of F is about (2−F)(1−F)pn with pn the frequency of the null allele. At low S, the expected bias on F̂ is therefore approximately twice, and that on Ŝ little less than four times, the frequency of the null allele (Figure 1). This bias can override the actual (null or small) F-value in outbred species. Note that similar biases are expected for other types of partial dominance problems (Supplementary Appendix 4). A general approximation for the bias on selfing rates is α(1−S)(2−S)/(1+α(1−S)), with α the fraction of heterozygotes read as homozygotes. The bias is of order 2α when S and α are low. For all sources of artifacts, including null alleles, the bias decreases when S increases.

Figure 1
figure 1

Difference, ΔS, between the selfing rate estimated from the inbreeding coefficient (F) and the actual selfing rate (S) as a function of S, when null alleles affect inbreeding estimates. The curves are based on equations provided in Supplementary Appendix 4. From top to bottom, the null-allele frequency is 0.1, 0.03 and 0.01, respectively.

Full size image

The influence of genotyping errors on PAA estimates has not been sorted out quantitatively, although Ritland notes in MLTR documentation that family estimates might be seriously biased by misscoring. We illustrate the case of null alleles using a simple model developed in Supplementary Appendix 6. We consider a specific (perhaps worst case) situation in which both null homozygous mothers and homozygous offspring are discarded from the analysis. In such a situation, the selfing rate is overestimated, and when S is small, the bias is about 4p/3 (where p is the frequency of the null allele). It decreases with increasing selfing rates (Figure 1 in Supplementary Appendix 6). Although this bias should not be neglected, it is much lower than that when S is estimated from the inbreeding coefficient. In addition, as mentioned above, null alleles result in inconsistencies between mother and offspring genotypes, and should therefore be more easily diagnosed in progeny arrays than in single-generation samples. The probability of detecting at least one critical offspring genotype (incompatible with mother genotype) rises fast with the number of offspring (n). For example, it can be shown from Table 1 in Appendix 6 that if the mother is heterozygous for two nonnull alleles this probability is about 2/3 when S is low and n=10. On the whole, genotyping errors more severely limit the inbreeding-coefficient approach than the PAA.

Getting estimates of S free from technical bias, especially under the PSA, therefore seems a worthwhile goal. This is exactly what is achieved by the multilocus method of David et al. (2007) that is insensitive to the technical artifacts that generate apparent heterozygote deficiencies, including null alleles. This is because the estimation of g2 is independent of F: technical biases create heterozygote deficiencies (bias in F), but do not create correlations in heterozygosity among loci because they occur independently at different loci. An illustration is provided in Figure 2. This method is particularly useful when S is low, that is when technical biases such as null alleles are more of a problem.

Figure 2
figure 2

Bias on single-generation estimates of S using the inbreeding coefficient method averaged over loci (top two curves) and the multilocus method based on g2 (bottom curves). Each value was obtained from 300 simulations of a sample of 100 individuals. Ten loci with gene diversity He=0.8 each were considered. Two misscoring rates were simulated, α=0.1 and α=0.3, with α the fraction of heterozygotes scored as homozygotes.

Full size image

Methodological bias

These biases are essentially related to violation of assumptions regarding population structure and history or to selection, and have variable consequences on Ŝ depending on the method (Table 1; Supplementary Appendix 6). Biases in Ŝ due to population structure and history are generally positive (Table 1). Various processes enter this category: (1) The first is biparental inbreeding, but its magnitude is likely to be limited. Jarne and Auld (2006) estimated from the plant data set from Goodwillie et al. (2005) that biparental inbreeding on average inflated apparent selfing rates by about 3%, irrespective of the true value. (2) The Wahlund effect results from mixing genetically differentiated subpopulations. Estimates of S are inflated because this process generates heterozygote deficiencies. If these populations differ in heterozygosity, the multilocus estimates obtained from the method of David et al. (2007) method will also be inflated. In general, the Wahlund effect is unlikely to play a significant role because sampling rarely includes strongly differentiated subpopulations. This is perhaps more questionable with the PAA, because allelic frequencies in the male gamete pool may differ from those in the female gamete pool (for example, as a consequence of pollinator behavior in plants). Such a difference can be evaluated using specific genetic analyses. (3) Deviation from inbreeding equilibrium is an issue for the PSA only, mainly when S is large and varied dramatically in the recent past. In the latter case the estimate will be between the current and past values of S. Such events can be detected using the multilocus method proposed by Enjalbert and David (2000). (4) Linkage disequilibrium is of course more of an issue for all multilocus methods, under both the PSA and the PAA, especially at high selfing rates. Under the PAA, it leads to overestimating S, and underestimating biparental inbreeding (Hedrick and Ritland, 1990). Its influence is less one-sided under the method of David et al. (2007) method. Whatsoever linkage equilibrium should be checked independently (see for example, Rousset and Raymond, 1997). The influence of asexual reproduction in those species able to both self-fertilize and reproduce clonally has not been sorted out for the multilocus PSA (David et al., 2007), and can formally be taken into account in the PAA. On the other hand, it might slightly lower estimates of the inbreeding coefficient (Prugnolle et al., 2005). The problem can partially be alleviated by discarding copies of mutilocus genotypes represented more than once.

Selection is a serious issue for all methods, as inbreeding depression is a prevalent feature of most species, including highly inbred ones (Husband and Schemske, 1996). The inbreeding coefficient at conception, corresponding to the primary selfing rate (value at fertilization before any selection), differs from that in adults as a consequence of inbreeding depression (see Supplementary Appendix 3). Whatever the method, estimates derived from adult individuals might therefore underestimate the primary selfing rate, especially in outcrossing species in which high inbreeding depression is expected early in ontogeny (Doums et al., 1996; Husband and Schemske, 1996). In this respect, it is worth conducting PSAs on individuals as young as possible. The same problem occurs in principle with the PAA, but is often mitigated in practice because progenies are genotyped as either seeds, or juveniles. Estimates can be corrected using independent estimates of inbreeding depression, accounting for selection up to the stages at which S is estimated. Another way out is the methodology proposed by Ritland (1990a): sampling at two stages (say seedlings and adults) allows jointly estimating F, S and inbreeding depression expressed between the two stages (Supplementary Appendix 3).

Computational problems can be encountered when the selfing rate is numerically estimated using maximum likelihood (multilocus PSA and PAA). Convergence is indeed not granted, or might be slow. Ritland (1990b, 2002) used two methods for maximizing likelihoods partly as a solution to this issue. On the other hand, estimating likelihoods allows building tests among competing models, a desirable quality for building a hypothesis-testing framework.

Comparing estimates derived from different methods

We have described several methods for quantifying inbreeding in populations, and one might be interested in comparing parameter estimates across methods. The main reason is that estimates from the PSA and the PAA represent different views of inbreeding, the first assesses inbreeding over several generations, while the second is more a ‘here and now’ estimate. Another more practical reason is that the maternal F is estimated together with S among offspring in the PAA.

We have already seen that estimates of S derived from the inbreeding-coefficient and multilocus (PSA) approaches might fruitfully be compared, for example to uncover potential technical pitfalls. This point is illustrated in David et al. (2007): these authors provide the striking example of an allozymic data set based on 12 populations of the freshwater snail Physa acuta. F-estimates suggested selfing rates equal or higher than 0.5, while the multilocus approach detected deviation from random mating in two populations only. This last result is much more in line with what is known from PAA and population structure analysis based on microsatellites (Henry et al., 2005). A more systematic comparison of results from the two approaches would be worthwhile, and might even provide insights on the magnitude of technical problems in species where the inbreeding rate is known from other approaches.

F-estimates of S have also been compared to PAA estimates (Jarne and Auld, 2006). It should be realized that a one-to-one relationship is not expected because estimates are not derived from the same life-history stage. Adult mothers should have lower apparent selfing rate than their offspring because of inbreeding depression (Supplementary Appendix 3). In the quite restricted animal data set, this difference was not observed, although inbreeding depression is usually strong in outbreeding species, suggesting that maternal F-estimates were overestimated due to technical problems (Jarne and Auld, 2006). The same relationship in plants is more in line with expectations, and consistent with some inbreeding depression at low selfing rates and more limited technical bias. However, one-third of F-estimates are negative, a surprising result with no obvious interpretation.

Recommendations

We propose simple recommendations for future studies (see also Table 1). These are mostly ‘rules of thumb’ based on the arguments developed above and our experience, providing estimates of selfing rates with a precision of, say, 10%. For more specific goals (for example, testing for difference among families or subpopulations), appropriate designs should be built.

(1) Markers: Preferentially use microsatellites, and several loci. A good target is LH=5, where L is the number of loci and H their average genetic diversity. (2) Sampling: n=30–50 individuals is a minimum for the PSA; larger numbers are required in the PAA, because sampling is structured in families (typical numbers are 20 families of 5–10 offspring, or more if maternal genotypes have to be inferred). There certainly is a trade off between the number of families and the number of offspring per family depending on whether the focus is on within-family or among-family variation. Individuals should be sampled as young as possible if one is interested in primary selfing rates. (3) Analysis: use both inbreeding coefficients and the multilocus method (RMES) when analyzing population samples, and multilocus maximum likelihood (MLTR) for the PAA. (4) Data interpretation should carefully consider how technical biases and violation of assumptions might interfere with the estimation of S. In this respect, several complementary analyses provide invaluable help, whatever the approach. These include evaluating genotyping repeatability, comparing single-locus and multilocus estimates, comparing different age-classes or estimating inbreeding depression, documenting the population genetic context, especially population substructure and gametic disequilibria, and jointly using the PAA and PSA whenever possible.

Conclusion and perspectives

Estimates of S can be derived from a variety of approaches, and it is certainly worth combining them because they document various aspects of the selfing rate, which is indeed a variable, dynamic parameter. No method is devoid of assumptions, and the interpretation of results from other analyses (for example, estimating inbreeding depression, analyzing population genetic substructure) has proved valuable in several species. Fortunately the scope of the PAA has been continuously enlarged over the years (Ritland, 2002), and the approach of David et al. (2007) renders the PSA as robust as the PAA. However, there is room for theoretical improvements, for example to account for partial dominance in the PAA (see Supplementary Appendix 6) or to extend its scope to other biological situations (see MLTR documentation). A significant effort should also be devoted to reducing technical biases as much as possible, as they arguably are the most important threat to the validity of selfing estimates.

Finally, the benefit of increasingly refining methods for estimating selfing is not to slowly converge toward the true selfing rate of a species, for there is no such thing. On the contrary, removing estimation bias and error should allow researchers to focus on what selfing really is: a variable trait and the product of a complex interaction between genotypes and environment (Brown et al., 1989), which constitutes a worthy subject of study on its own. Even flowers on the same plant can display different selfing rates, and the selfing rate may also vary during the lifetime of an individual (for example, along the reproductive season in a snail). Properly quantifying fine-grained variation in S is a prerequisite to understanding how much of this variation reflects adaptation versus environmental or developmental constraints. A promising avenue for further methodological developments is to refine the decomposition of the variance in selfing rates among different levels of biological organization (populations, families, individuals or even different ramets or successive reproductive events of a single individual) and their covariance with relevant environmental variables in natural as well as in experimental settings.