Rare sex or out of reach equilibrium? The dynamics of FIS in partially clonal organisms

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{F_{IS}}<0 $$\end{document}FIS¯<0. Yet in partially clonal species, both FIS < 0 and FIS > 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 state-and-time discrete Markov chain model to describe the dynamics of FIS 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 FIS = 0. Thus, although clonality alone does not lead to departures from Hardy-Weinberg expectations once reached the final equilibrium state, both negative and positive FIS 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 FIS and reduces its rate of change over time, leading to a hyperbolic increase of the maximal time needed to reach the final mean \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{F_{IS,\infty }} $$\end{document}FIS,∞¯ value expected at equilibrium. Conclusion Our results argue for a dynamical interpretation of FIS 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{F_{IS}} $$\end{document}FIS¯ values deviating from zero and/or a large variation of FIS over loci are observed. Electronic supplementary material The online version of this article (doi:10.1186/s12863-016-0388-z) contains supplementary material, which is available to authorized users.


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][8][9][10][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.
The interpretation of standard population genetic indices from partially clonal populations can be challenging, as expectations are likely to depend on the rate of clonality, which is usually unknown in natural populations. The estimate of this rate on the basis of indirect approaches such as population genetics analysis remains elusive. One example of an index that has been suggested to change with the rate of clonality is F IS [13,14]. Within diploid individuals, it represents a correlation coefficient among alleles at a particular locus, and depends on their tendency to be randomly associated F IS = 0 or more likely identical (F IS > 0) or not identical (F IS < 0) at the population level. F IS is defined either based on population heterozygosity (H eexpected heterozygosity, H oobserved heterozygosity) or allele identities/homozygosity (Fallele identity within individuals, Θallele identity within the population [13]): 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 Hardy-Weinberg 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 Hardy-Weinberg 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 F IS;∞ À ¼ −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 dynamically-stable 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], 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 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.

Mathematical model
The biological template for our model is a single population with a finite number of individuals. These individuals correspond to ramets, i.e. factually or potentially physiologically distinct units that may or may not be genetically identical or descended from the same parent. All individuals follow the same life cycle, which consists of a dominant diploid phase during which they can acquire heritable mutations (Fig. 1). All individuals are monoecious and can produce offspring both by clonal and sexual (here defined as random mating including selfing) reproduction. A fixed number of these offspring individuals corresponding to a constant population size survive randomly to replace their parents in the next generation. 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): I Mutation. The theoretical frequencies ν ii,I , ν ij,I of each genotype after mutation are derived as: where α = 1 − μ, the probability that an allele does not mutate, and β ¼ μ n−1 , the probability that an allele mutates into one of the n − 1 others during one generation. This corresponds to a classic k-alleles or Jukes-Cantor substitution model [25]. II Gamete formation (allele segregation then fusion). The gamete frequencies in the gamete pool after sexual reproduction ν i,I are calculated as: There is no difference in the allele frequencies between sexes, mating types etc., and all individuals contribute equally to the gamete pool (pangamy).
III Reproduction. The genotype frequencies ν ii,III , ν ij,IIIafter reproduction are calculated as: based on the results from Eqs. 1 and 2. The rate of clonality c thus corresponds to the proportion of offspring per generation that is the result of clonal reproduction. The remainder of the offspring ("rate of sexuality", 1 − c) is derived from random mating including autogamy, assuming panmixis.
8 < : Fig. 1 Schematic overview of the mathematical model (example for two alleles). In a dominantly diploid population of fixed size N, the number of individuals/ramets q with a certain genotype (here aa, aA, or AA) at a particular locus, observed at generation t, and the corresponding genotype frequencies ν = q/N may change due to mutation (here symmetrical from a to A and from A to a with rate μ; see Eq. 1), reproduction (random mating at rate 1 − c; see Eqs. 2 and 3) or genetic drift (modeled by multinomial drawing of N individuals from the genotype frequency distribution; see Eq. 4), until observation at the next generation IV Genetic Drift. The vector of genotype frequencies v tþ1 → at the next generation depending on the population size is derived from: where X t+1 is the state of the model at the next generation, drawn from a multinomial distribution M that is based on N, the population size counting all potentially reproducing individuals (mathematically the number of samples), and ν → III , the vector of genotype frequencies derived from Eq. 3 (mathematically the probabilities of the genotype "categories"). Transition probabilities P between any two model states X t , X t + 1 can then be calculated based on: where q ii, t + 1 , q ij, t + 1 ∈ ℕ 0 are the counts (natural numbers) of individuals per genotype in the presumed next state X t+1 which sum to N. Note that our description of genetic drift is based on genotype frequencies rather than allele frequencies. As explained in [22], describing population genetic processes based on allele frequencies is a mathematical convenience justified by HWE (i.e. assuming exclusively sexual reproduction), which assures that allele frequencies can always be directly translated into genotype frequencies.
For partially clonal populations, we cannot automatically assume HWE and thus modeled all population genetic processes, including genetic drift, at the genotypic level.

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 q → tþ1 ¼ 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.
Secondly, we approached the full model by following the dynamics of F IS over time from three different start states for combinations of c, μ and N representative of the different regions of the parameter space. Aggregating the transition probabilities between all model states in a transition matrix M (same current state per column, i.e. columns summing to one), the probability distribution of the model states (and consequently the probability distribution of F IS ) at time t, given by the vector x → t , is derived by matrix multiplication: where x → 0 describes the start state (vector of zeros except for a single one at X 0 ). The steady state distribution of population states corresponds to the dominant eigenvector of the transition matrix M. Based on the transition matrix M, we also calculated the time to the steady state distribution (Markov chain mixing time, Additional file 1, 1.5).
We illustrated the numerical result of our model using three start states: , 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][29][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 F IS ;∞ e (Markov chain mixing time, Additional file 1, 1.5).
To link our results with those obtained by previous authors, we calculated the final mean F IS ;∞ À from Eq. 10 in [13] formalized for metapopulation: , q s is the probabilities that two individuals taken at random in the same reproductive subpopulation after migration were sired in the same reproductive subpopulation one generation before, and q d the probability that two individuals taken at random in different reproductive subpopulations after migration originated from the same subpopulation one generation before. By setting q s = 1, q d = 0 (a single finite population, no migration) and s = 1/N (random mating): We derived the corresponding expected convergence time of F IS À iteratively obtained from Eq. 5 in [13] (detailed in Additional file 1, 1.6) 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-/homo-zygosity as continuous variables whatever the population size considered.
Finally, to get a better idea how our theoretical results are comparable to those published for field data, we looked at the sampling effect of using different numbers of polymorphic loci L to estimate the mean F IS À of the population at time t, F IS;t;L À . Assuming that each locus represents an independent estimate of this mean (no confounding effect of linkage), and that the genotype frequencies are known exactly (exhaustive sampling of all individuals/ramets), it is derived as: Both assumptions are usually violated [14,18], so that our results represent a conservative estimate of the true error of this method. We randomly sampled both the steady state distribution F IS;∞ e and the instantaneous distribution F IS ;50 e of a population that was 50 generations ago at HWE at all loci with equal allele frequencies (isoplethy for two alleles per locus), for the same parameter combinations that we previously used to illustrate the dynamics of the full model. Based on 10 5 random samples of size L, we then calculated the mean signed deviation of the sample means F IS,t ⋅ from the true mean F IS;t À : where z 1 and z 2 represent the number of positive and negative deviations, respectively. Loci at or near fixation are typically not used in population genetic studies, since they are especially affected by genotyping errors. We therefore excluded all loci where the frequency of one allele exceeds 1− ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1= 2N ð Þ p (near fixation; Additional file 2: Figure S2.1 for the derivation of this value, and compare similar considerations in [31]) from the calculation of values for this analysis.
All computations were performed in Python 2.7 with 64 bit precision, using the modules numpy, scipy [32], networkx [33] and matplotlib [34] (basic code in Additional file 3). We illustrate some of our results with de Finetti diagrams [24,35] (Fig. 2a, c, e), which are ternary plots of the genotype frequencies [ν aa , ν aA , ν AA ] at one locus with two alleles within a population (see Additional file 2: Figures S2.1 to S2.4 for more information). Details for the literature review in the discussion are given in Additional file 4.

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 (self-incompatibility 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 maximal convergence time due to reproduction only, t c (derived in Additional file 1, 1.2) can be approximated as with ε ¼ 1 2N corresponding to a small error term as convergence to the steady state distribution (HWE: Fig. 2a, Additional file 2: Figure S2.2) is asymptotic. Between t c = 1 for exclusively sexual and t c = ∞ for exclusively clonal populations, the dependence of t c on c is not linear, but shapes like an hyperbola of type y ¼ 1 þ constant x (Fig. 2b). Consequently, the time The maximal convergence time due to mutation only t μ (derived in Additional file 1, 1.2) can be approximated as which simplifies for two alleles to 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 μ ¼ 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.
The maximal convergence time due to reach a steady state distribution (Fig. 2e, Additional file 2: Figure S2.4) due to genetic drift only t N depends on the population size and on the number of initial genotypes (Fig. 2f ). It grows approximately linearly with N and the slope of this dependence σ increases with the number of genotypes/alleles (Additional file 1, 1.2): with τ ≥ 1 and 1 < σ < 2 determined from numerical results (Additional file 2: Table S2.2) and simulations (Additional file 2: Table S2.3). For example, in a population of 100 individuals, up to about 160 generations can be required until only one genotype remains at a biallelic locus, and loci with more alleles take even a few generations more. The lowest maximal convergence time amongst t c , t μ , t N , can serve as approximation to determine the evolutionary process that dominates the dynamics of genotype frequencies and its associated steady state. Pairwise equality between the convergence times can thus be sued to partition the parameter space of our model distribution (eqs. 1 to 4, c ∈ [0, 1], μ∈]0, 0.5], N ∈ [1, ∞ [) into three domains where either process dominates (Fig. 3). Based on eqs. 12 to 14, these pairwise equalities resolve to: 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.
Keeping population size and mutation rate constant (N = 100, μ = 10 − 6 ) while successively increasing the rate of clonality illustrates the changes in the dynamics of F IS as genetic drift takes over (Fig. 4a to e, Additional file 2: Figures S2.6 and S2.7, Table S2.1). At low to intermediate rates of clonality (here 0 < c ≤ 0.8), the dynamics of F IS appear almost identical to those expected for a purely sexual population ( Fig. 4a and b). However, variation around the final mean  Fig. 4d). Indeed, transition probabilities that decrease F IS raise while those that increase F IS reduce (compare transition probabilities from (25,50,25) to (25,49,26) and to (25,51,24) in Additional file 2: Table S2.1). Moreover, probabilities to stay on the same genotypic state increase with clonality, except for states near fixation when mutation is high. 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. 4c-e). In parallel, the probability of allele fixation becomes increasingly dependent on the intial F IS,0 value ( Fig. 4a-e, Additional file 2: Figures S2.6 and S2.7, left columns). The first reason for this effect, which also explains the convergence to negative F IS;∞ À 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. 4a-e). 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 F IS;∞ À ¼ −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 F IS;t (starting from F IS,0 = 0 in sexual populations) has not reached the final F IS;∞ À 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:  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) Population size N = 100 throughout. Columns: crate of clonality, F IS,∞ = 0mutation rate, t Ngenetic drift maximal expected convergence time, t creproduction maximal convergence time, t μmutation maximal convergence time, t Iconvergence time to the mean F IS;∞ À based on the model in [13], t IIconvergence time to the mean F IS;∞ À based on our model, t IIIconvergence time to full final distribution of F IS;∞ e . Rows: example parameter sets (compare Fig. 4). Bold: min (t c , t μ ) Thus, when comparing different rates of clonality in populations with the same size and realistic mutation rate, the increase of the convergence time to F IS;∞ À is directly related to t c until t c > t μ .
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) For example in highly clonal populations (Fig. 4d, e at t = 50), the average deviations from the current exact mean 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 dimensiontimeto the description of the steady state distribution F IS ;∞ e and its mean F IS;∞ À 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 Hardy-Weinberg 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 F IS;∞ À ¼ 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 F IS;∞ À 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 long-lasting 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 F IS;t;L À unless a high number of loci is used. These findings thus contribute to explaining the frequent reports of negative mean 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 F IS;∞ À , 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 F IS À in natural populations of partially clonal organisms.

The respective effects of reproduction, mutation and genetic drift in partially clonal populations
Results on the dynamics of genotype frequencies due to each evolutionary process (model parameters) formally demonstrate that: 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 F IS;∞ À 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). F IS;∞ À ¼ 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 F IS;∞ À 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. F IS;∞ À ¼ 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 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 F IS;∞ À 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).
The deceleration of the dynamics of F IS during its approach to the steady state is connected to the hyperbolical increase of t c (Fig. 2b) and thus much stronger under high rates of clonality. We demonstrated that a comparison of the maximal expected convergence times t c , t μ , t N can be an efficient means to predict the overall pattern of F IS dynamics (Figs. 3, 4 and 5, Additional file 2: Figure S2.5). The times t c and t μ can even be used to estimate convergence times of the complete model ( can be estimated by the minimum of t c and t μ (i.e. usually t c in small populations). If t N ≪ min(t c , t μ ), loci with different initial genotype frequencies may not converge to the same final F IS,∞ value (convergence to genotypic uniformity), so that the expected final F IS;∞ À ¼ −1 cannot be reached within biologically realistic time spans even under nearly pure clonality (Fig. 4e, Additional file 2: Figure S2.5).
Though not yet included in our model, perenniality leading to overlapping generations (partial survival of the individuals across generations) is expected to slow down F IS dynamics even further. If disturbances are sufficiently frequent, e.g. in very instable environments or in populations cyclically changing between exclusive sexual and clonal reproduction [43,44], the final F IS;∞ À and F IS;∞ e based on the currently observed rate of clonality may even never be reached.
Finally, during the convergence toward the final F IS ;∞ e , 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 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 F IS ;∞ e , 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 non-exhaustive sampling, genotyping errors (e.g. undetected null alleles for SSRs) or preferential sampling of loci with near-isoplethy (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).
We kept studies on species with self-incompatibility 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 F IS ;∞ e . 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 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 F IS;∞ À under equilibrium as the reference. For example, F IS;t;L À ¼ −0:083 the exhaustively sampled population 16 ( Fig. 5)  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 F IS;t;L À over loci, the full distribution of F IS ;t;L e 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
Our results allow reconciling predictions for F IS under partial clonality from theoretical models, which suggest departures from F IS = 0 only at nearly pure clonality, with some empirically observed values, which show such departures also where sex is known or suspected to be frequent. Examining the dynamic effects of clonality under the varying influence of mutation and drift showed three main implications for interpreting F IS in partially clonal populations: non-negative 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 DNA-based 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) e F IS;t , exact distribution of F IS at a time t F IS;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, Hardy-Weinberg 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