Rare sex or out of reach equilibrium? The dynamics of F _{ IS } in partially clonal organisms
 Katja Reichel^{1},
 JeanPierre Masson^{1},
 Florent Malrieu^{2},
 Sophie ArnaudHaond^{3} and
 Solenn Stoeckel^{1}Email author
https://doi.org/10.1186/s128630160388z
© The Author(s). 2016
Received: 29 January 2016
Accepted: 1 June 2016
Published: 10 June 2016
Abstract
Background
Partially clonal organisms are very common in nature, yet the influence of partial asexuality on the temporal dynamics of genetic diversity remains poorly understood. Mathematical models accounting for clonality predict deviations only for extremely rare sex and only towards mean inbreeding coefficient \( \overline{F_{IS}}<0 \). Yet in partially clonal species, both F _{ IS } < 0 and F _{ IS } > 0 are frequently observed also in populations where there is evidence for a significant amount of sexual reproduction. Here, we studied the joint effects of partial clonality, mutation and genetic drift with a stateandtime discrete Markov chain model to describe the dynamics of F _{ IS } over time under increasing rates of clonality.
Results
Results of the mathematical model and simulations show that partial clonality slows down the asymptotic convergence to F _{ IS } = 0. Thus, although clonality alone does not lead to departures from HardyWeinberg expectations once reached the final equilibrium state, both negative and positive F _{ IS } values can arise transiently even at intermediate rates of clonality. More importantly, such “transient” departures from Hardy Weinberg proportions may last long as clonality tunes up the temporal variation of F _{ IS } and reduces its rate of change over time, leading to a hyperbolic increase of the maximal time needed to reach the final mean \( \overline{F_{IS,\infty }} \) value expected at equilibrium.
Conclusion
Our results argue for a dynamical interpretation of F _{ IS } in clonal populations. Negative values cannot be interpreted as unequivocal evidence for extremely scarce sex but also as intermediate rates of clonality in finite populations. Complementary observations (e.g. frequency distribution of multiloci genotypes, population history) or time series data may help to discriminate between different possible conclusions on the extent of clonality when mean \( \overline{F_{IS}} \) values deviating from zero and/or a large variation of F _{ IS } over loci are observed.
Keywords
Background
Reproductive systems impact the evolution of genetic diversity at the population level [1, 2], making them an important factor for considerations on the evolvability of species. Partially clonal species, i.e. species that are able to reproduce both sexually and clonally, are common across many phyla and ecosystems [3] and represent an important part of the global biodiversity. They include many species of which evolutions are directly impacted by, or impacting humans, such as cultivated species [4], pathogens [5], invasive species [6], and species threatened by extinction (e.g. [7–11]). Partially clonal species are therefore frequently the subject of molecular analyses describing their genetic diversity [12], and the conclusions drawn depend on a correct understanding of the effects of their reproductive mode on the genetic composition of their populations.
Results from both definitions differ only for loci with just a single allele remaining (fixation), where F _{ IS } cannot be defined.
To date, only few mathematical models studying F _{ IS } at selectively neutral loci in partially clonal populations have been published. For partially clonal populations otherwise complying with the HardyWeinberg conditions, F _{ IS } and the underlying genotype frequencies are expected to be identical to those expected for random mating, except for the approach to the HardyWeinberg equilibrium (HWE) being slowed down as the rate of clonality increases [15]. If mutation and genetic drift are taken into account [13], very high rates of sexual reproduction are expected to lead eventually to strongly negative mean F _{ IS } values up to \( \overline{F_{IS,\infty }}=1 \) for completely clonal populations. In addition to this effect on the mean, the shape of the full steady state (i.e. “equilibrium”: the convergence of values toward a given probabilistic and dynamicallystable distribution) distribution of F _{ IS }, measured by its variance, skewness and kurtosis, also changes with the rate of clonality [16]. Based on the results of [13], \( \overline{F_{IS}} \) was suggested as an informative parameter to estimate the rate of clonality [14, 17] in connection with other indices such as linkage disequilibrium or the frequency of repeated multiloci genotypes [18]. However, using the mean of the steady state distribution provided by [13] as a reference for the mean \( \overline{F_{IS}} \) values from field studies often pointed to rates of clonality that were at odds with other indices or even direct observation (e.g. [19, 20]).
While some previous works may be interpreted as demonstrating negative F _{ IS } as a signature of nearly exclusive clonality [13, 14], others underline the influence of clonality not only on the steady state distributions of F _{ IS } but also on the temporal dynamics of this parameter in natural population [15, 16]. Here we aimed at complementing the results of the previous studies by describing the temporal changes of genotype frequencies over time under the distinct and joint influence of partial clonality, mutation and genetic drift. Understanding how quickly the steady state distribution of F _{ IS } is reached again after a disturbance (e.g. change of reproductive system, stochastic or selective events) due to clonality, mutation and drift, may indeed help to explain the unexpected departures observed in nature, in populations otherwise considered as undergoing rather frequent sexual reproduction.
We used a stochastic model to follow the neutral dynamics of genotype frequncies in the basic case of a single locus in a diploid, isolated and panmictic population that combines random mating and clonality. We derived the dynamical effects of population size, mutation rate and rate of clonality separately on genotype frequencies over generations and how long it takes until a steady state distribution is reached again after any disturbance. We subsequently connected these partial results to analyze the “complete” system with the joint effects of reproductive mode, mutation and genetic drift. Finally, we discuss how our results may provide new hypotheses in the interpretation of field data, based on examples from a literature review, and we provide methodological recommendations for future analyses of genotype frequencies in partially clonal populations.
Methods
Mathematical model
We translated this system into a time and state discrete Markov chain model, conceptually similar to [16] to follow the population dynamics of genotype frequencies at one locus with multiple alleles, as in the classic models of Wright [21] and others [22, 23]. Each time step of the model corresponds to one generation, i.e. the time between two consecutive observations of the population (Fig. 1). The model states represent all possible distributions of the N individuals on g genotypes (a genotypic state): For a single locus n alleles, the number of genotypes corresponds to g = n(n + 1)/2 and the number of states to (N + g − 1) !/(N ! ⋅ (g − 1) !). For example, a locus of 3 alleles (A1, A2 and A3) results in 6 possible genotypes (A1A1, A1A2, A1A3, A2A2, A2A3, A3A3) and a population composed of 100 individuals could evolve among the 96 560 646 single possible repartitions of 100 individuals within the 6 genotypes (hereafter named a population state) [16, 24].
At each time step, the population makes a transition from its current state to a next state (where current and next state can be the same), based on a vector of transition probabilities. These probabilities depend on the genotype frequencies ν _{ ii }, ν _{ ij } (with the indices i ≠ j denoting different alleles) derived from the current population state, and on the three constant model parameters: population size N, mutation rate μ and rate of clonality c, according to the following equations (Fig. 1; simplified for the special case of two alleles in Additional file 1, 1.1):
There is no difference in the allele frequencies between sexes, mating types etc., and all individuals contribute equally to the gamete pool (pangamy).
Model analysis and description of biological consequences
First, we studied the effect of each of the three model parameters (c, μ, N) on the genotype frequencies by itself. Setting the other two parameters to have no influence on the model result, i.e. c = 1, μ = 0 and/or N = ∞ (or no random drawing in Eq. 4), and substituting Eq. 2 into Eq. 3, the model reduces to one equation per process, i.e. Eq. 1 for μ, Eq. 3 (with 2) for c, and Eq. 4 for N. For each equation/process, we then determined the steady state distributions of genotype frequencies at one locus, i.e. the probabilistic combination of population states for which \( {\overrightarrow{q}}_{t+1}={\overrightarrow{q}}_t \). We also derived the respective maximal convergence times t _{ c }, t _{ μ } and t _{ N } of each separate evolutionary process to reach steady state distributions of genotype frequencies. While t _{ N } could only be approximated from numerical results (Markov chain first passage time approach and simulations), for c and μ convergence to the steady states is asymptotical as it can be described by geometric progressions (details of derivation in Additional file 1, 1.2). We therefore defined a universal “acceptable error” ɛ = 1/(2N), corresponding to one half the minimal change in genotype frequency that would be measurable by exhaustive sampling in a population of finite size N, below which the distance from the steady states has to pass (convergence criterion). Using the maximal convergence times t _{ c }, t _{ μ } and t _{ N } as measures for the “strength” with which each process acts upon the genotype frequencies, we could then use this analytical basis to partition the parameter space of the full model into regions where either process dominates the genotype frequency dynamics.
We illustrated the numerical result of our model using three start states: F _{ IS,0} ∊ { − 1; 0; 1} under isoplethic allele frequencies (equally frequent, \( {\nu}_a={\nu}_A=\frac{1}{n} \)), standing for HWE (F _{ IS,0} = 0) and the most extreme deviations from it (complete homozygosity, F _{ IS,0} = 1; complete heterozygosity, F _{ IS,0} = − 1). These start states were chosen to represent the range of F _{ IS } values, and not because of their biological meaning or how they may be reached in nature. Deviations from steady state distributions may derive from a recent change in the rate of clonality (e.g. from full sexuality with F _{ IS,0} = 0), or full adaptation to past selection for (F _{ IS,0} = − 1) or against (F _{ IS,0} = 1) heterozygotes, or changes in population size (demographic bottleneck, founder event) [26, 27], or secondary contact between two populations in which different alleles become fixed (F _{ IS,0} = 1) or hybridization with subsequent reproductive isolation from the parents (F _{ IS,0} = − 1) (e.g. [28–30]). Based on the transition matrix M, we also calculated the time to the final distribution of states (= steady states), which is also the time until the exact final distribution of \( \tilde{{F_{IS}}_{,\infty }} \) (Markov chain mixing time, Additional file 1, 1.5).
Where F and θ are the allele identities within individuals and between individuals respectively. In contrast to our own model, these equations do not contain the number of alleles. This is because they are based on an infinite alleles model, and treat the expected and observed hetero/homozygosity as continuous variables whatever the population size considered.
Examples for empirical F _{ IS } values of partially clonal populations from published field studies
To get a summary view of the kind of F _{ IS } values that can be encountered in real field populations suspected to reproduce using partial clonality, we compiled the data for 51 populations belonging to 13 species (seven angiosperms, four protists, a sponge and a nematode), based on 13 previously published studies (see references in Additional file 3) selected for their near fit with the assumptions of our model from a Web of Science search for [(microsatellite OR “SSR” OR “simple sequence repeat” OR “SNP” OR “single nucleotide polymorphism”) AND (clonal OR asexual OR vegetative OR apomictic OR apomixis OR agamospermy OR parthenogenesis)]. All studies are based on SSR data. When F _{ IS } values were not directly provided by publications, we calculated F _{ IS } values per locus from the reported H _{ o } and H _{ e }. We grouped populations into three classes according to the information given by the authors about their putative rate of clonality: i) rarely clonal, ii)frequent clonality and sexuality (including unknown) and iii)rarely sexual. We included populations for which preferential outbreeding (selfincompatibility system) was expected.
Results
The maximal times until the genotypic diversity at one locus has converged to its steady state distribution changed greatly with the rates of clonality. During this convergence, the steady state distributions, the dynamics of F _{ IS } and the underlying genotype frequencies through time were also affected by the rate of clonality. However, both maximal convergence times and the dynamics of genotypic diversity at one locus strongly depend on the interactions of the rate of clonality with mutation rate and population size. Indeed, each of these evolutionary processes has its own maximal convergence time toward the steady state distribution of genotype frequencies and F _{ IS } (t _{ c }, t _{ μ }, t _{ N }). The overall convergence time of a population to its steady state distribution, which is shaped by mutation, genetic drift and the reproduction mode, depends on the relative “strength” of each of those three forces.
The steady state distribution for mutation corresponds to F _{ IS } = 0 and, due to the symmetry of mutation between alleles in our models, equal allele frequencies (Additional file 1, 1.3; Additional file 2: Figure S2.3). For \( \mu =\frac{1}{n} \), in theory the highest mutation rate (each allele has the same chance to mutate or not to mutate into other alleles), t _{ μ } is only one generation. For realistic mutation rates ranging from 10^{−3} to 10^{−18} [36, 37], convergence due to mutation can take much longer: setting ɛ = 0.005 and assuming a mutation rate of μ = 10^{−6}, t _{ μ } would be around 2.6 ⋅ μ ^{−1} or 2.6 million generations.
If t _{ c } ≪ t _{ μ }, t _{ N }, as is usually the case for strictly sexually reproducing populations, genotype frequencies quickly converge to HWE, i.e. F _{ IS,∞} ≅ 0. In very large and highly clonal populations where , genotype frequencies also converge, although more slowly, to HWE so that t _{ μ } ≪ t _{ c } , t _{ N } is equally awaited in this case. This contrasts with small and highly clonal populations where F _{ IS,∞} ≅ 0, as dominant genetic drift does not lead to convergence toward t _{ N } ≪ t _{ μ } , t _{ c }, but rather to a successive loss of genotypes and eventually to genotypic uniformity (fixation of one homozygous or heterozygous genotype F _{ IS,∞} = 0 or − 1). In conclusion, both random mating and mutation lead to a random association of alleles (Additional file 2: Table S2.1, [16]), so that the dynamics of F _{ IS } in (partially) clonal populations is mostly driven by the relative “strength” of genetic drift.
A closer look at the transitions between the predominance of either process shows that changes are actually gradual. This is because different processes don’t globally compensate each other, as each convergence pattern is different (Fig. 2a, c, e, Additional file 2: Table S2.1). We therefore examined the joint action of all three processes in more details.
For loci starting with a high genotypic diversity, a marked difference in the positive and negative ranges of the F _{ IS } distribution, namely the progressive disappearance of positive values, appears after some generations (t ~ 50, Fig. 4ce). In parallel, the probability of allele fixation becomes increasingly dependent on the intial F _{ IS,0} value (Fig. 4ae, Additional file 2: Figures S2.6 and S2.7, left columns). The first reason for this effect, which also explains the convergence to negative \( \overline{F_{IS,\infty }} \) in highly clonal populations, is the trend towards randomly increasing the frequency of the heterozygous genotype(s) at HWE and approximately evenly distributed allele frequencies (Additional file 1, 1.4). Indeed, HWE is a concave function in the F _{ IS } space (for one locus, 2 alleles, see Additional file 2: Figure S2.1) which in probability facilitates slightly negative F _{ IS } as a result of stochastic changes even for random mating populations (see Additional file 2: Figure S2.1 and Table S2.1, e.g. transition probabilities from (25,50,25) to (25,49,26) and (25,51,24)). Partial and full clonality only prevents populations, partly or completely, to get back to HWE, thus to express this feature over more generations. The second reason is the dependence of the rate of fixation (of one allele) on the initial F _{ IS,0} at extreme clonal rate (Fig. 4ae). All loci reaching positive F _{ IS,t } values then have a higher chance to move quickly toward fixation, helped by each random mating event (as the HWE parabola in de Finetti diagrams connects fixation states) and thus to be usually excluded from F _{ IS } calculation. Over time, this only leaves analyzable loci with negative values. It can be noted that contrastingly, loci shifting to negative F _{ IS,t } (higher probabilities than shifting to positive, due to the concave shape of HWE) are unlikely to reach \( \overline{F_{IS,\infty }}=1.0 \) by pure stochastic genetic drift, as each sexual and homoplasy events bring back populations toward the HWE concave parabola (Fig. 2a, c, e). Interestingly in our example, the mean \( \overline{F_{IS,t}} \) (starting from F _{ IS,0} = 0 in sexual populations) has not reached the final \( \overline{F_{IS,\infty }} \) even after 200 generations while t _{ N } = 159 because convergence to drift steady states is slowed down by “weak” mutation.
Highly clonal populations dominated by mutation rather than genetic drift present a different picture (Fig. 4f, g: c = 1.0, μ = {10^{− 2}, 10^{− 1}}, N = 100, or simulation for c = 1.0, μ = 10^{− 3}, N = 10^{3} in Additional file 2: Figure S2.5). As in predominantly sexually reproducing populations, F _{ IS } values converge to only slightly negative final \( \overline{F_{IS,\infty }} \) (Fig. 4g: \( \overline{F_{IS,\infty }}\approx 0.02 \)) but at a speed that mainly depends on t _{ μ } and with a limited variation of F _{ IS } values across generations. In contrast, the instantaneous \( \overline{F_{IS,t}} \) distributions appears more symmetrical as the fixation of single alleles is very rare (Fig. 4f, g). As the maximal allelic diversity (number of possible alleles) increases beyond two, t _{ μ } increases accordingly and mutation is “weakened” even further in comparison to the other processes (Fig. 2d).
Convergence times under partial clonality for seven representative example parameter sets
N = 100  c  μ  t _{ N }  t _{ c }  t _{ μ }  t _{ I }  t _{ II }  t _{ III } 

A  0.0  10^{−6}  159  1  2.6 × 10^{6}  1  1  2.6 × 10^{6} 
B  0.8  10^{−6}  159  25  2.6 × 10^{6}  27  27  2.6 × 10^{6} 
C  0.97  10^{−6}  159  159  2.6 × 10^{6}  174  177  2.6 × 10^{6} 
D  0.99  10^{−6}  159  529  2.6 × 10^{6}  464  498  2.6 × 10^{6} 
E  1.0  10^{−6}  159  ∞  2.6 × 10 ^{ 6 }  38,366  ≫ 40,000  2.6 × 10^{6} 
F  1.0  10^{−2}  159  ∞  263  234  138  264 
G  1.0  10^{−1}  159  ∞  25  25  14  25 
Finally, an important observation is also that when approaching the steady state in partially clonal populations, the variations of F _{ IS } increase again as compared to variations expected at the steady states, as the “evolutionary memory” of clonality (sensu [38]) ease the influence of genetic drift on genotype frequencies (Fig. 4). Partial clonal populations that have not yet reached their steady state distribution of \( \tilde{{F_{IS}}_{,\infty }} \) may thus show greater deviations from their current exact mean \( \overline{F_{IS,t}} \) value than expected for the final \( \overline{F_{IS,\infty }} \) values, especially in small population sizes (Additional file 2: Figures S2.8 and S2.9). For example in highly clonal populations (Fig. 4d, e at t = 50), the average deviations from the current exact mean \( \overline{F_{IS,50}} \) based on ten loci even exceeded ±0.1 the variation yet expected at the steady states.
Discussion
Our results on the dynamics of F _{ IS } in partially clonal populations add a new dimension – time – to the description of the steady state distribution \( \tilde{{F_{IS}}_{,\infty }} \) and its mean \( \overline{F_{IS,\infty }} \) derived from previous models [13, 16]. They allow estimating the time needed to reach steady state distributions of genotype frequencies, be there are different or not from the classical HardyWeinberg proportions defined for sexual populations. Those estimates not only take into account the hyperbolic relationship between the rate of clonality and the convergence times to HWE [15], but also the effects of two other major evolutionary forces, genetic drift and mutation (Fig. 3).
Thus far, the frequent occurrence of departure from HWE observed in partially clonal species was interpreted as fitting a narrow range of scenario where clonality (heterozygote excesses) and/or selfing (homozygote excesses) would dominate. This study offer a new perspective on these observations, by associating significant F _{ IS } to a much broader range of stable and transient scenario including modest rates of clonality, even in the absence of selfing. These results thus call for increased caution when interpreting these field data alone to make inferences on the importance of clonality or selfing, in the absence of other biological knowledge of the mating system.
Indeed, our work confirms the previous findings of departure from the usual average \( \overline{F_{IS,\infty }}=0 \) only at very high rates of clonality [13, 14, 17]. However, it also recalls and supports the pioneering work from Marshall and Weir [15] on the speed of convergence towards this \( \overline{F_{IS,\infty }} \) in partially clonal populations. We extend the previous results by elucidating how the interplay between rate of clonality, population size (i.e. genetic drift) and mutation affects the nature of the steady state, the time required to reach it and the transient variation of the distribution of F _{ IS } before reaching it. The results reported here show that owing to an increased “evolutionary memory” (sensu [38]) of past genotypic diversity in partially clonal populations, population history such as changes in reproductive mode can produce a transient but possibly longlasting overrepresentation of F _{ IS } values departing from 0, even under modest rates of clonality. Additionally the variation of F _{ IS }, which changes dynamically until the steady state is reached, increases the risk to misestimate the mean \( \overline{F_{IS,t,L}} \) unless a high number of loci is used. These findings thus contribute to explaining the frequent reports of negative mean \( \overline{F_{IS}} \) in the literature even for species expected to undergo clonal reproduction at intermediate rates. Their otherwise restricted interpretation as a signature of extreme clonality, conceived under the narrow light of the steady state distributions of \( \overline{F_{IS,\infty }} \), is thus alleviated. In the following, we further discuss the way this increased “evolutionary memory” and the interaction between rate of clonality, genetic drift and mutation affect the interpretation of \( \overline{F_{IS}} \) in natural populations of partially clonal organisms.
The respective effects of reproduction, mutation and genetic drift in partially clonal populations

Compared to the standard expectation for exclusively sexual populations, clonality only slows down the approach to HWE (Fig. 2b), at a rate depending on other processes such as mutation and genetic drift to which it grants an increased influence.

mutation, if acting independently at each allelic copy, leads towards F _{ IS,∞} = 0 (Fig. 2c and Additional file 2: Figure S2.3) even for extreme clonal rates, when genetic drift is comparatively negligible (μ > 1/N)

genetic drift, if not negligible, tends towards genotypic uniformity, i.e. either. F _{ IS,∞} = − 1 (one heterozygous genotype remains) or fixation (one homozygous genotype remains), depending on the initial genotype frequencies (Fig. 2e).
Thus the negative \( \overline{F_{IS,\infty }} \) predicted by previous models, and observed in many empirical studies, strongly depend not only on the rate of clonality but also on its interplay with the population size and, to a lesser extent, the mutation rate (Fig. 3). \( \overline{F_{IS,\infty }}=0 \) is expected at steady state even under high rates of clonality and pure clonality in large populations (Additional file 2, 2.6), as well as in smaller populations when mutation rates are high, as e.g. for microsatellites markers (Fig. 3, “mutation” part of the diagram; Fig. 4f, g). Contrastingly, negative \( \overline{F_{IS,\infty }} \) as suggested by [13] are also expected, yet only in smaller clonal populations or when back mutations are neglected (as in [13]), thus only where genetic drift dominates the dynamics of genotype frequencies (Fig. 4d, e; Fig. 3, “genetic drift” part of the diagram).
The maximal convergence time for mutation, t _{ μ }, to its steady state situated on the HWE (i.e. \( \overline{F_{IS,\infty }}=0 \)) increases with the number of possible alleles (Fig. 2d). This steady state due to mutation exists regardless of the number of alleles or of the mutation scheme as long as all alleles mutate at constant rates (e.g. [39, 40], rates need not be equal among alleles, compare Additional file 1, 1.3). This result is only violated if mutation depends on the other alleles at the same locus within the same individual, as for gene conversion [41, 42]. In this particular case, which we did not assess here but which promotes homozygote excess, populations should converge with a similar dynamics to F _{ IS } = 1 or to the fixation of one allele.
According to our results, clonal populations dominated by mutation and with no gene conversion (Fig. 3, “mutation” part of the diagram) may distinguish themselves from their counterparts dominated by random mating (Fig. 3, “reproduction” part of the diagram) by the rarity of loci that are fixed for one allele throughout the whole population, rather than by a different mean \( \overline{F_{IS}} \). In order to assess the relative importance of clonality, the examination of the distribution of F _{ IS } values among multiple loci and of the proportion of fixed loci would thus be more informative than the mean F _{ IS } value itself.
The dynamics of F _{ IS } in partially clonal populations: the implications of the long lasting temporal dynamics and of the large interloci variance
Our equations and numerical results demonstrated that the dynamics of genotype frequencies and F _{ IS } are slowed down in partially clonal populations, which therefore retain traces of their past for much longer than their exclusively sexual counterparts. This questions the generic value of the final mean \( \overline{F_{IS,\infty }} \) as derived in previous studies [13, 16] as a basis for the interpretation of field data. For at least intermediate rates of clonality, the genetic and genotype composition of the population may indeed reflect population history rather than the present day reproductive system and thus mostly depend on the time since the last disturbance of the population (Fig. 4, Additional file 2, 2.6).
Finally, during the convergence toward the final \( \tilde{{F_{IS}}_{,\infty }} \), loci may go through intermediate distributions that would be definitely unusual in exclusively sexual populations both in terms of genotype frequencies and F _{ IS } (Fig. 4, Additional file 2: Figure S2.6). Even at modest rates of clonality, the variation of F _{ IS } is increased compared to exclusively sexual populations. Consequently, information from more loci is required to reach an accurate estimate of the mean \( \overline{F_{IS,t}} \) in partially clonal populations compared to strictly sexual ones. We found that the variation of F _{ IS } observed during the approach to the steady state distribution may be even greater than predicted based on the final \( \tilde{{F_{IS}}_{,\infty }} \), especially when t _{ N } ≫ t _{ μ } (Additional file 2: Figures S2.8 and S2.9).
Transient F _{ IS } values: a new hypothesis to account for values observed in field data?
We performed a literature review to illustrate the frequent observation of a very wide variety of F _{ IS } values, positive as well as negative, in partially clonal populations (Fig. 5; details in Additional file 4). Field data may be influenced by technical biases, including sampling bias due to an unknown spatial structure of clones, missing rare genotypes due to nonexhaustive sampling, genotyping errors (e.g. undetected null alleles for SSRs) or preferential sampling of loci with nearisoplethy (thus increasing the probability to find negative F _{ IS }). Moreover, biological processes other than those in our model may have acted on the data we collected. We therefore applied strict criteria to standardize the dataset we used and retain only those populations that fitted our model best. We present only studies that did include repeated multiloci genotypes in their calculations, published F _{ IS } values (or H _{ o } and H _{ e }) per locus and clearly isolated population, and reported data on organisms whose life cycle fits with our model (i.e. dominantly diploid, no cyclic clonality as e.g. in aphids). Only few studies matched these strict criteria [20, 45–56]. We kept studies on species with selfincompatibility systems, as this system of preferential outbreeding has been shown to have very little effect on F _{ IS } at loci physically (nearly) unlinked to the SI genes [57, 58].
Owing to the increased evolutionary memory of past demographic fluctuations, the results of our study open up new possible explanations for the presence, but also the absence, of positive and negative F _{ IS } values (both at individual loci or the mean) in partially clonal populations. On their way back to the steady state distribution, even intermediately clonal populations may transiently exhibit F _{ IS,t } values that should be rare in the final steady state distribution \( \tilde{{F_{IS}}_{,\infty }} \). Such values can be due to the increased variation of F _{ IS } or echo the past departure from equilibrium due to population history (e.g. demographic bottleneck, change in the rate of clonal reproduction). As an example of how to apply our results, slightly negative mean \( \overline{F_{IS,t,L}} \) over loci in some wild cherry populations (populations 14–16 in Fig. 5, [20]), would have suggested almost exclusive clonality when taking the expected mean \( \overline{F_{IS,\infty }} \) under equilibrium as the reference. For example, \( \overline{F_{IS,t,L}}=0.083 \) the exhaustively sampled population 16 (Fig. 5) with N = 247, μ ∈ [10^{− 3}, 10^{− 12}] suggest c ≈ 0.98. However, the proportion of repeated multiloci genotypes and inferences from parentage analysis suggested an intermediate rate of clonality instead (c ~ 0.5) [20, 59]; based on this value, genotype dynamics would be dominated by random mating (t _{ c } ≈ 10). Using our theoretical results, the observed F _{ IS } distribution and its negative mean could be explained either by extensive logging in the past (population history) or by the fact that only nine loci were analyzed in small populations of about a hundred individuals (increased variation; compare Additional file 2, 2.8 and 2.9).
Our results also suggest ways to further improve the population genetic inferences in natural populations of partial clonals, based on F _{ IS } in connection with other parameters as proposed in [14]. First, other parameters than F _{ IS } may be much more informative when attempting to assess the relative importance of sexual versus clonal reproduction in the dynamics and evolution of partially clonal population. Maximizing the number of loci studied by moving from population genetics to population genomics may help to improve the statistical basis of inferences of population parameters. Second, if pursuing the investigation of the influence of rates of clonality on F _{ IS }, rather than focusing exclusively on the mean \( \overline{F_{IS,t,L}} \) over loci, the full distribution of \( \tilde{{F_{IS}}_{,t,L}} \)values per population should be reported and interpreted. Collecting time series of samples may also provide valuable information, as field data normally represent only a “snapshot” of genotype frequencies at a particular point in time, that may or may not be representative of the steady state distribution of the parameters being studied. Using the Markov chain model implemented here, it is not only possible to statistically analyze example trajectories, but also to analytically derive the transition probabilities between two consecutive sets of genotype frequencies for a range of clonal rates, based on population size, mutation rate and number of generations between time series samples. The results of this study open up perspectives for the development of a unified statistical method to infer rates of clonal reproduction, or other population genetic parameters of our model if c would be known, based on the analysis of temporal samples. Taking into account the temporal dynamics of genetic descriptors of the populations, including F _{ IS } would therefore help to improve the biological interpretation of values from field data, or to refine methods for estimating the rate of clonality based on a collection of population genetic indices.
Conclusions

nonnegative F _{ IS }, including null values, are not a reliable indicator of the absence of clonal reproduction, as they may occur i) even at steady state under highly frequent but not exclusive clonality, provided the influence of mutation compared to drift in populations is large (μ ≫ 1/N). Or they may occur ii) transiently under all rates of clonality, a likely situation for many wild populations considering the hyperbolic relationship between the rate of clonality and the time toward convergence reported here.

negative F _{ IS } values are not a reliable indicator of nearly exclusive clonal reproduction, as significant deviations from F _{ IS } = 0 for multiple loci are also expected after departures from the steady state distribution in partially clonal populations, which generally last longer even if the rate of clonality is only intermediate.

An increased number of loci are required to maintain the accuracy of DNAbased estimates of population genetic parameters in partially clonal compared to exclusively sexual populations in general, and the study of time series rather than single snapshots of genetic data may lead to more accurate estimates of the rate of clonal reproduction in particular.
Nomenclature
F _{ IS }, inbreeding coefficient represents a correlation coefficient among alleles at a particular locus within polyploid (diploid here) individuals
t, current generation (discrete time)
\( \overset{\sim }{{\mathrm{F}}_{\mathrm{IS},\mathrm{t}}} \), exact distribution of F _{ IS } at a time t
\( \overline{{\mathrm{F}}_{\mathrm{IS},\mathrm{t}}} \), mean F _{ IS } at a time t
H_{o}, observed heterozygosity
H_{e}, expected heterozygosity
F, allele identity within individuals
Θ, allele identity within the population
HWE, HardyWeinberg equilibrium
N, population size
n, number of alleles at one locus. Alleles can be named as A, a for a biallelic locus or A1, A2, …Ai for locus with more than two alleles
i ≠ j ≠ k ≠ l, indices referring to alleles
v _{ i } , allele frequency of Ai
g, number of different genotypes at one locus
q_{ij}, number of individual of genotype AiAj
v _{ij}, genotype frequency of genotype AiAj
c, rate of clonal reproduction
s, rate of selfing, set here at 1/N as expected under random mating
μ, mutation rate
α, the probability that an allele does not mutate
β, the probability that an allele mutates into one of the n − 1 others during one generation
X, random variable
P(X), probability of a random variable X
M, transition matrix
ℳ, multinomial distribution
L, numbers of studied polymorphic loci
t_{c}, maximal number of generations to convergence due to rate of clonality
t_{μ}, maximal number of generations to convergence due to mutation rate
t_{N}, maximal number of generations to convergence due to genetic drift
min(t_{c}, t_{μ}), minimum convergence time between two evolutionary forces
ε, universal “acceptable error” corresponding to one half the minimal change in genotype frequency that would be measurable by exhaustive sampling in a population of finite size N
q_{s}, probabilities that two individuals taken at random in the same reproductive subpopulation after migration were sired in the same reproductive subpopulation one generation before
q_{d}, probability that two individuals taken at random in different reproductive subpopulations after migration originated from the same subpopulation one generation before
Declarations
Acknowledgements
We thank Jurgen Angst, Valentin Bahier, Fabien Halkett, Cédric Midoux, Stéphane De Mita, Martin Lascoux and two anonymous reviewers for helpful discussions and comments on the manuscript. We greatly profited from the exchange of ideas within the CLONIX group.
Authors thank the French National Research Agency, CLONIX project (ANR11BSV70007), the French National Institute for Agricultural Research, Plant Health and the Environment Department (INRASPE) and the Région Bretagne, France for supporting this research.
Funding
This work was supported by the French National Research Agency, CLONIX project (ANR11BSV70007). Katja Reichel worked in this project during her PhD granted by the French National Institute for Agricultural Research, Plant Health and the Environment Department (INRASPE) and the Région Bretagne, France.
Availability of data and materials
The data supporting the results of this article are included within the article and its Additional files. A basic version of the model code is available as Additional file 3 and from http://www6.rennes.inra.fr/igepp_eng/Productions/Software.
Authors’ contributions
JPM and SS laid the foundation for this work. JPM, FM, KR and SS conceived the mathematical model. KR and SS performed the scientific programming and generated the numerical data. KR did the literature review, most of the data analyses and all illustrations. KR, SAH and SS wrote the manuscript. All authors contributed to editing, managed and performed by SS. SAH and SS were responsible for funding and grant applications. All authors have read and approved the final manuscript.
Competing interests
The authors declare that they have no competing interests.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Duminil J, Fineschi S, Hampe A, Jordano P, Salvini D, Vendramin GG, Petit RJ. Can population genetic structure be predicted from lifehistory traits? Am Nat. 2007;169:662–72.View ArticlePubMedGoogle Scholar
 Duminil J, Hardy OJ, Petit RJ. Plant traits correlated with generation time directly affect inbreeding depression and mating system and indirectly genetic structure. BMC Evol Biol. 2009;9:177.View ArticlePubMedPubMed CentralGoogle Scholar
 de Meeûs T, Prugnolle F, Agnew P. Asexual reproduction: Genetics and evolutionary aspects. Cell Mol Life Sci. 2007;64:1355–72.View ArticlePubMedGoogle Scholar
 McKey D, Elias M, Pujol B, Duputié A. The evolutionary ecology of clonally propagated domesticated plants: Tansley review. New Phytol. 2010;186:318–32.View ArticlePubMedGoogle Scholar
 Tibayrenc M, Ayala FJ. Reproductive clonality of pathogens: A perspective on pathogenic viruses, bacteria, fungi, and parasitic protozoa. Proc Natl Acad Sci. 2012;109:E3305–13.View ArticlePubMedPubMed CentralGoogle Scholar
 Liu J, Dong M, Miao SL, Li ZY, Song MH, Wang RQ. Invasive alien plants in China: role of clonality and geographical origin. Biol Invasions. 2006;8:1461–70.View ArticleGoogle Scholar
 Luijten SH, Oostermeijer JGB, van Leeuwen NC, den Nijs HC. Reproductive success and clonal genetic structure of the rare Arnica montana (Compositae) in the Netherlands. Plant Syst Evol. 1996;201:15–30.View ArticleGoogle Scholar
 Sydes MA, Peakall R. Extensive clonality in the endangered shrub Haloragodendron lucasii (Haloragaceae) revealed by allozymes and RAPDs. Mol Ecol. 1998;7:87–93.View ArticleGoogle Scholar
 Brzosko E, Wróblewska A, Ratkiewicz M. Spatial genetic structure and clonal diversity of island populations of lady’s slipper (Cypripedium calceolus) from the Biebrza National Park (northeast Poland). Mol Ecol. 2002;11:2499–509.View ArticlePubMedGoogle Scholar
 Setsuko S, Ishida K, Tomaru N. Size distribution and genetic structure in relation to clonal growth within a population of Magnolia tomentosa Thunb. (Magnoliaceae). Mol Ecol. 2004;13:2645–53.View ArticlePubMedGoogle Scholar
 Brzyski JR, Culley TM. Genetic variation and clonal structure of the rare, riparian shrub Spiraea virginiana (Rosaceae). Conserv Genet. 2011;12:1323–32.View ArticleGoogle Scholar
 Schön I, Martens K, Dijk P, editors. Lost Sex. Dordrecht: Springer Netherlands; 2009.Google Scholar
 Balloux F, Lehmann L, de Meeûs T. The population genetics of clonal and partially clonal diploids. Genetics. 2003;164:1635–44.PubMedPubMed CentralGoogle Scholar
 Halkett F, Simon J, Balloux F. Tackling the population genetics of clonal and partially clonal organisms. Trends Ecol Evol. 2005;20:194–201.View ArticlePubMedGoogle Scholar
 Marshall DR, Weir BS. Maintenance of genetic variation in apomictic plant populations. Heredity. 1979;42:159–72.View ArticleGoogle Scholar
 Stoeckel S, Masson JP. The exact distributions of F _{ IS } under partial asexuality in small finite populations with mutation. PLoS One. 2014;9:e85228.View ArticlePubMedPubMed CentralGoogle Scholar
 de Meeûs T, Lehmann L, Balloux F. Molecular epidemiology of clonal diploids: a quick overview and a short DIY (do it yourself) notice. Infect Genet Evol. 2006;6:163–70.View ArticlePubMedGoogle Scholar
 ArnaudHaond S, Duarte CM, Alberto F, Serrão EA. Standardizing methods to address clonality in population studies. Mol Ecol. 2007;16:5115–39.View ArticlePubMedGoogle Scholar
 Motoie G, Ferreira GEM, Cupolillo E, Canavez F, PereiraChioccola VL. Spatial distribution and population genetics of Leishmania infantum genotypes in São Paulo State, Brazil, employing multilocus microsatellite typing directly in dog infected tissues. Infect Genet Evol. 2013;18:48–59.View ArticlePubMedGoogle Scholar
 Stoeckel S, Grange J, FernándezManjarres JF, Bilger I, FrascariaLacoste N, Mariette S. Heterozygote excess in a selfincompatible and partially clonal forest tree species  Prunus avium L. Mol Ecol. 2006;15:2109–18.View ArticlePubMedGoogle Scholar
 Wright S. Evolution in Mendelian populations. Genetics. 1931;16:97–159.PubMedPubMed CentralGoogle Scholar
 Ewens WJ. Mathematical population genetics: I. Theoretical Introduction. 2nd ed. New York: Springer; 2004.View ArticleGoogle Scholar
 Gale JS. Theoretical population genetics. London. Springer: Unwin Hyman; 1990.View ArticleGoogle Scholar
 Reichel K, Bahier V, Midoux C, Masson JP, Stoeckel S. Interpretation and approximation tools for big, dense Markov chain transition matrices in ecology and evolution. Algorithms Mol Biol. 2016;10:1–14.Google Scholar
 Jukes TH, Cantor CR. Evolution of protein molecules. In: Munro HN, editor. Mammalian protein metabolism. Volume III. New York: Academic; 1969. p. 21–132.View ArticleGoogle Scholar
 Nei M. Molecular population genetics and evolution. Amsterdam: North Holland Publishing Company; 1975.Google Scholar
 Brakefield PM. The variance in genetic diversity among subpopulations is more sensitive to founder effects and bottlenecks than is the mean: a case study. In: Fontdevila A, editor. Evolutionary biology of transient unstable populations. Berlin: Springer Berlin Heidelberg; 1989. p. 145–61.View ArticleGoogle Scholar
 Edmands S, Feaman HV, Harrison JS, Timmerman CC. Genetic consequences of many generations of hybridization between divergent copepod populations. J Hered. 2005;96(2):114–23.View ArticlePubMedGoogle Scholar
 Zalapa JE, Brunet J, Guries RP. The extent of hybridization and its impact on the genetic diversity and population structure of an invasive tree, Ulmus pumila (Ulmaceae). Evol Appl. 2010;3(2):157–68.View ArticlePubMedPubMed CentralGoogle Scholar
 Lindtke D, Buerkle CA, Barbará T, Heinze B, Castiglione S, Bartha D, Lexer C. Recombinant hybrids retain heterozygosity at many loci: new insights into the genomics of reproductive isolation in Populus. Mol Ecol. 2012;21:5042–58.View ArticlePubMedGoogle Scholar
 Graffelman J, Camarena JM. Graphical tests for HardyWeinberg equilibrium based on the ternary plot. Hum Hered. 2008;65:77–84.View ArticlePubMedGoogle Scholar
 Oliphant TE. Python for scientific computing. Comput Sci Eng. 2007;9:10–20.View ArticleGoogle Scholar
 Hagberg AA, Schult DA, Swart PJ. Exploring network structure, dynamics, and function using NetworkX. In Proceedings of the 7th Python in Science Conference (SciPy2008), SciPy Edited by: Varoquaux G, Vaught T, Millman J. CA USA: Pasadena; 2008:11–15.Google Scholar
 Hunter JD. Matplotlib: a 2D graphics environment. Comput Sci Eng. 2007;9:90–5.View ArticleGoogle Scholar
 de Finetti B. Considerazioni matematiche sull’ereditarietà mendeliana. Metron. 1926;6:3–41.Google Scholar
 Drake JW, Charlesworth B, Charlesworth D, Crow JF. Rates of spontaneous mutation. Genetics. 1998;148:1667–86.PubMedPubMed CentralGoogle Scholar
 Hile SE, Yan G, Eckert KA. Somatic mutation rates and specificities at TC/AG and GT/CA microsatellite sequences in nontumorigenic human lymphoblastoid cells. Cancer Res. 2000;60:1698–703.PubMedGoogle Scholar
 Desai MM. Reverse evolution and evolutionary memory. Nat Genet. 2009;41:142–4.View ArticlePubMedGoogle Scholar
 Estoup A, Jarne P, Cornuet JM. Homoplasy and mutation model at microsatellite loci and their consequences for population genetics analysis. Mol Ecol. 2002;11:1591–604.View ArticlePubMedGoogle Scholar
 Ellegren H. Microsatellites: simple sequences with complex evolution. Nat Rev Genet. 2004;5:435–45.View ArticlePubMedGoogle Scholar
 McMahill MS, Sham CW, Bishop DK. Synthesisdependent strand annealing in meiosis. PLoS Biol. 2007;5:e299.View ArticlePubMedPubMed CentralGoogle Scholar
 Flot JF, Hespeels B, Li X, Noel B, Arkhipova I, Danchin EGJ, Hejnol A, Henrissat B, Koszul R, Aury JM, Barbe V, Barthélémy RM, Bast J, Bazykin GA, Chabrol O, Couloux A, Da Rocha M, Da Silva C, Gladyshev E, Gouret P, Hallatschek O, HecoxLea B, Labadie K, Lejeune B, Piskurek O, Poulain J, Rodriguez F, Ryan JF, Vakhrusheva OA, Wajnberg E, et al. Genomic evidence for ameiotic evolution in the bdelloid rotifer Adineta vaga. Nature. 2013;500:453–7.View ArticlePubMedGoogle Scholar
 Allen DE, Lynch M. The effect of variable frequency of sexual reproduction on the genetic structure of natural populations of cyclical parthenogen. Evolution. 2012;66:919–26.View ArticlePubMedGoogle Scholar
 Berg LM, Lascoux M. Neutral genetic differentiation in an island model with cyclical parthenogenesis. J Evol Biol. 2000;13:488–94.View ArticleGoogle Scholar
 Duran S, Pascual M, Estoup A, Turon X. Strong population structure in the marine sponge Crambe crambe (Poecilosclerida) as revealed by microsatellite markers. Mol Ecol. 2004;13:511–22.View ArticlePubMedGoogle Scholar
 Nagamitsu T, Ogawa M, Ishida K, Tanouchi H. Clonal diversity, genetic structure, and mode of recruitment in a Prunus ssiori population established after volcanic eruptions. Plant Ecol. 2004;174:1–10.View ArticleGoogle Scholar
 Corral JM, Molins MP, Aliyu OM, Sharbel TF. Isolation and characterization of microsatellite loci from apomictic Hypericum perforatum (Hypericaceae). Am J Bot. 2011;98:e167–9.View ArticlePubMedGoogle Scholar
 Jiang K, Gao H, Xu NN, Tsang EPK, Chen XY. A set of microsatellite primers for Zostera japonica (Zosteraceae). Am J Bot. 2011;98:e236–8.View ArticlePubMedGoogle Scholar
 Tew JM, Lance SL, Jones KL, Fehlberg SD. Microsatellite development for an endangered riparian inhabitant, Lilaeopsis schaffneriana subsp. recurva (Apiaceae). Am J Bot. 2012;99:e164–6.View ArticlePubMedGoogle Scholar
 Liu W, Zhou Y, Liao H, Zhao Y, Song Z. Microsatellite primers in Carex moorcroftii (Cyperaceae), a dominant species of the steppe on the QinghaiTibetan Plateau. Am J Bot. 2011;98:e382–4.View ArticlePubMedGoogle Scholar
 Barnabe C, Buitrago R, Bremond P, Aliaga C, Salas R, Vidaurre P, Herrera C, Cerqueira F, Bosseno MF, Waleckx E, Breniere SF. Putative panmixia in restricted populations of Trypanosoma cruzi isolated from wild Triatoma infestans in Bolivia. PLoS One. 2013;8:e82269.View ArticlePubMedPubMed CentralGoogle Scholar
 Tesson SVM, Borra M, Kooistra WHCF, Procaccini G. Microsatellite primers in the planktonic diatom Pseudonitzschia multistriata (Bacillariophyceae). Am J Bot. 2011;98:e33–5.View ArticlePubMedGoogle Scholar
 Villate L, Esmenjaud D, Van Helden M, Stoeckel S, Plantard O. Genetic signature of amphimixis allows for the detection and fine scale localization of sexual reproduction events in a mainly parthenogenetic nematode. Mol Ecol. 2010;19:856–73.View ArticlePubMedGoogle Scholar
 Gao H, Jiang K, Geng Y, Chen XY. Development of microsatellite primers of the largest seagrass, Enhalus acoroides (Hydrocharitaceae). Am J Bot. 2012;99:e99–101.View ArticlePubMedGoogle Scholar
 McInnes LM, Dargantes AP, Ryan UM, Reid SA. Microsatellite typing and population structuring of Trypanosoma evansi in Mindanao, Philippines. Vet Parasitol. 2012;187:129–39.View ArticlePubMedGoogle Scholar
 Vilas R, Cao A, Pardo BG, Fernández S, Villalba A, Martínez P. Very low microsatellite polymorphism and large heterozygote deficits suggest founder effects and cryptic structure in the parasite Perkinsus olseni. Infect Genet Evol. 2011;11:904–11.View ArticlePubMedGoogle Scholar
 Navascués M, Stoeckel S, Mariette S. Genetic diversity and fitness in small populations of partially asexual, selfincompatible plants. Heredity. 2009;104:482–92.View ArticlePubMedPubMed CentralGoogle Scholar
 Stoeckel S, Klein EK, OddouMuratorio S, Musch B, Mariette S. Microevolution of Sallele frequencies in wild cherry populations: respective impacts of negative frequency dependent selection and genetic drift. Evolution. 2012;66:486–504.View ArticlePubMedGoogle Scholar
 Stoeckel S. Impact de la propagation asexuée et du système d’autoincompatibilité gamétophytique sur la structuration et l’évolution de la diversité génétique d’une essence forestière entomophile et disséminée, Prunus avium L. France: Cemagref; 2006.Google Scholar