## Abstract

Genome duplication (polyploidy) is a recurrent evolutionary process in plants, often conferring instant reproductive isolation and thus potentially leading to speciation. Outcome of the process is often seen in the field as different cytotypes co-occur in many plant populations. Failure of meiotic reduction during gametogenesis is widely acknowledged to be the main mode of polyploid formation. To get insight into its role in the dynamics of polyploidy generation under natural conditions, and coexistence of several ploidy levels, we developed a general gametic model for diploid–polyploid systems. This model predicts equilibrium ploidy frequencies as functions of several parameters, namely the unreduced gamete proportions and fertilities of higher ploidy plants. We used data on field ploidy frequencies for 39 presumably autopolyploid plant species/populations to infer numerical values of the model parameters (either analytically or using an optimization procedure). With the exception of a few species, the model fit was very high. The estimated proportions of unreduced gametes (median of 0.0089) matched published estimates well. Our results imply that conditions for cytotype coexistence in natural populations are likely to be less restrictive than previously assumed. In addition, rather simple models show sufficiently rich behaviour to explain the prevalence of polyploids among flowering plants.

## 1. Introduction

Polyploidy (i.e. the multiplication of entire chromosome sets above the diploid level) has played a key role in evolution and diversification of many eukaryotic organisms, including yeasts, insects, vertebrates and plants [1]. Genome multiplication is especially common in monilophytes (horsetails and ferns) and angiosperms, with estimated proportions of species with one or more polyploidization events in their evolutionary history of up to 95 per cent and 70 per cent, respectively [2]. Such events can be relatively recent (i.e. neopolyploidization), but genomic data have shown that virtually all angiosperms underwent one or more rounds of ancient whole genome duplication [3]. Because polyploidization, unlike other speciation processes, may confer instant reproductive isolation, it is considered a major mode of sympatric speciation, at least in the plant kingdom [4]. Quantitative estimates of the contribution of polyploidy to speciation have only recently begun to appear; it has been estimated that 15 per cent of speciation events in angiosperms and 31 per cent in ferns are associated with ploidy increases [5]. Considering the ubiquity and significance of polyploidy, understanding this mode of speciation basically means understanding the single most common evolutionary transition in the evolutionary history of angiosperms.

Whereas polyploid derivates may originate through somatic chromosome doubling, unreduced gamete formation is undoubtedly the major driving force in the formation of polyploids [6]. The origin via the production of unreduced gametes has been originally applied to autopolyploids, but it may also be valid for allopolyploids that have relatively frequently been hypothesized to originate as diploid homoploid hybrids that polyploidize afterwards. The frequency of recent autopolyploids in natural conditions is co-determined by the rates of their formation, establishment success, longer-term persistence, reproductive success and further spread. In addition to studies that have addressed one or more of these factors experimentally [7], the last few decades have seen several theoretical models aiming to evaluate the fates of new polyploids within the populations of their diploid/lower-ploid ancestors [8–15]. Most models are based on the assumption that the newly formed polyploid must cope with the minority disadvantage (i.e. the frequency-dependent process in which the rare ploidy is evolutionary disadvantaged unless it develops breeding barriers or shows higher fitness or fertility [16]). However, some data indicate that odd-numbered interploidy hybrids (such as triploids) may actually enhance the rate of polyploid formation and facilitate rather than diminish their fixation [10]. The outcome of these two opposing, but not mutually exclusive, processes (i.e. triploid block versus triploid bridge) is largely affected by the reproductive and fitness characteristics of triploid crosses. In general, the fate of a newly arisen tetraploid in a diploid population depends on the rate of triploid formation (or, conversely, the strength of triploid block), triploid fitness (viability and fertility) relative to diploids and tetraploids, and the ploidy of functional gametes of triploids. Simulations have shown that even partially fit triploids can increase the likelihood of diploid–tetraploid coexistence and, in some cases, facilitate tetraploid fixation [10].

Despite the development of more complex models that take into account some key parameters that may affect the formation and establishment of polyploids, including ploidy fitness, ploidy fertility and the rate of unreduced gamete production, their outcomes were not checked against actual ploidy frequencies observed under natural conditions. Consequently, the relative role of triploid block versus triploid bridge cannot be assessed for existing field systems. Analysis of actual ploidy frequencies was for a long time impractical because the number of cytotyped samples in conventional karyological studies was limited by methodological constraints that precluded large-scale population screening. Broad insights into ploidy variation *in situ* only became possible with the advent of flow cytometry, which transformed plant karyology from the examination of single individuals to the analysis of whole populations across multiple species and different spatial scales [17]. Ploidy data for several thousand individuals are now commonplace in cytogeographic studies, which provide an ideal framework for relating *in silico* results of ploidy dynamics to empirical population data.

In this paper, we develop a simple and general gametic model of ploidy-level frequencies, on the basis of the models of Husband [10] and Yamauchi *et al.* [12]. Such models predict equilibrium frequencies of individual ploidy levels using a few basic parameters (probabilities of formation of unreduced gametes and fertilities of higher ploidy levels). In our approach, we have developed different versions of the model for systems comprising different ranges of ploidy levels. We first present an analysis of the model's behaviour for systems comprising individuals ranging from: (i) 2*x* to 3*x*, (ii) 2*x* to 4*x*, and (iii) 2*x* to 6*x*, exploring the model's general predictions in light of current knowledge of the dynamics of polyploid series and other models. We then fit the model to observed ploidy frequency data to: (i) find out whether the empirical frequencies are compatible with such a simple gametic generating process, and (ii) obtain estimates of these key parameters of the model and to compare their numerical values with independent estimates available from other sources.

Although there are two basic modes of polyploidy establishment (i.e. auto- and allopolyploidy), for the sake of simplicity, we deal only with autopolyploidy (also called intraspecific or intrapopulation polyploidy) to avoid problems related to the involvement of different parental species in the polyploid's origin. Recent data show that autopolyploids are more common than previously considered and have important evolutionary and ecological roles in natural populations [18].

## 2. Material and methods

### (a) A gametic model of polyploidy: general framework

Let us assume that the ploidy level of an individual fully determines its gametic structure. We further assume that diploids produce mostly reduced gametes and a certain proportion of unreduced gametes (*b*); triploids produce pairs of reduced (*x* and 2*x*) gametes, in addition to a certain proportion of unreduced gametes (*c*). We assume that aneuploid gametes either do not form or that aneuploid progeny does not survive.

We further assume that populations of individual cytotypes are fully homogeneous (i.e. gametic structure of all individuals of a given cytotype is identical) and that gametes combine freely in individuals of the same ploidy level (irrespective of the actual mode of reproduction, whether outcrossing or selfing). Thus, we do not assume any *a priori* minority cytotype disadvantage. The progeny of each individual are thus fully determined by the probabilities of individual gametic types (see [12]). (We refer to this later as the ‘gametic process’.) For example, a diploid individual that produces haploid gametes with probability 1-*b* and unreduced gametes with probability *b*, will have diploid progeny if two reduced gametes combine (i.e. with the probability (1 − *b*)^{2}), triploid with the probability 2*b*(1 − *b*), and tetraploid with the probability *b*^{2}. Using the same logic, the progeny of a triploid plant with *c* proportion of unreduced gametes will range from diploid (if two reduced gametes with the haploid chromosome number merge, with probability (1 − *c*)^{2}/4) to hexaploid (if two fully unreduced gametes merge, with probability *c*^{2}; see also electronic supplementary material, table S2 and figure S2).

We further assume that individuals of all ploidy levels reproduce only generatively, have identical lifespans and survival, and have non-overlapping generations. For the sake of simplicity, we also disregard the fact that hybridization between different ploidy levels can take place. In addition, we assume that the system is closed and that the maximum viable ploidy level is 6*x* (plant groups with ploidy variation ranging from diploid to above hexaploid are quite rare under natural conditions). Changes in numbers of individuals of ploidy *i* from time *t* to time *t* + 1 can then be approximated using a set of recurrent equations:
2.1where *P _{i}*(

*t*) is the abundance of plants of ploidy

*i*(the subscript thus equalling the particular multiple of the basic chromosomal complement, and in the general case ranging from 2 to 6) at the time

*t*,

*F*is the proportion of progeny with ploidy

_{ij}*i*of mother individuals of ploidy

*j*as determined by the gametic process,

*μ*is the relative fertility of individuals of

_{i}*i*th ploidy; and

*s*and

*n*are, respectively, survival to maturity of and number of seeds produced by a diploid plant. For convenience, we assume, without loss of generality,

*μ*

_{2}= 1.

*F _{ij}*'s follow directly from the gametic process outlined above (see also the electronic supplementary material). For example,
2.2a
2.2b
2.2c
2.2dwhere

*b*is the proportion of unreduced gametes in the diploid and

*c*is the proportion of unreduced gametes in the triploid.

The equilibrium abundances of individual ploidies are then given by the dominant eigenvector of the matrix of coefficients of the recurrent equation (2.1): 2.3

which depends only on the values of *μ _{i}* and

*F*. The latter in turn depends on probabilities associated with the formation of unreduced gametes in diploids (

_{ij}*b*) and triploids (

*c*).

### (b) Ploidy-level frequency data

Ploidy estimates using DNA flow cytometry were performed as a part of other projects during 2000–2011. In addition, data on species and population ploidy composition were excerpted from published articles. In order to guarantee comparability of the results, only species meeting the following criteria were included in the study: (i) the occurrence of at least two ploidy levels from the range between 2*x* and 6*x*, (ii) ploidy data available for a total of at least 90 individuals from at least three populations, (iii) the existence of a sexual mode of reproduction, and (iv) (presumed) autopolyploid origin of higher ploidy levels. Most data originated from contact zones or were collected across the entire distributional range of a species. In addition to species-level data, individual mixed-ploidy populations within these taxa with more than 90 cytotyped individuals were analysed separately, and are referred to here as *species_name* P#. While population-level data provide more accurate estimates of ploidy frequencies, they can suffer from random shifts in ploidy composition owing to stochastic processes. Our total dataset consisted of 26 angiosperm species and additional 13 populations within these species, and included ploidy data for nearly 40 000 individuals (see the electronic supplementary material, table S1). Species lacking diploids and having all ploidy levels divisible by one small integer (*Cerastium* spp., *Elytrigia* spp.) were analysed so that the lowest ploidy level was assumed to be diploid.

### (c) Fitting the model to the data: choice of models to be fit

We fit the gametic model to empirical data for each species/population. For the fitting, we introduced some simplifying assumptions. First, we assumed that only diploids and triploids are capable of producing unreduced gametes. Further, unless otherwise stated, we assumed that all *μ*_{3} = *μ*_{4} = *μ*_{5} … , i.e. fertility of all higher ploidy to be equal over all odd and even ploidy levels. Although fertility in triploids and pentaploids is often markedly lower [7,19], analysis of the model behaviour shows that the effect of *μ*_{3} and *μ*_{5} on model predictions is rather low relative to the fertility and survival of even-ploid plants.

For species with ploidy levels up to tetraploid, we fit a model with the parameter *c* (i.e. proportion of unreduced gametes in a triploid) set to zero (‘model I’, which is essentially the model of Husband [10]; table 1). The model thus has only two parameters, and the 2*x*–3*x*–4*x* system has two degrees of freedom, meaning that the system is just-identified. Further, to examine whether the assumption that *c* = 0 is reasonable in a system that contains all three ploidy levels, we relaxed the condition *c* = 0, re-fit the model and checked the value of *c*. For species having ploidy levels only up to triploid, we also fit model I with the additional constraint, *μ*_{3} = 0 (referred to as the ‘minimal model’). In species/populations with ploidy levels from 2*x* up to 6*x*, we fit the model with parameters *b, c* and *μ* (‘model II’).

### (d) Fitting the model to the data: techniques

Parameters of model I and the minimal model were determined exactly by finding their analytic solutions by hand or using Mathematica v. 8 (Wolfram Research, Inc., Champaign, IL, USA; see the electronic supplementary material S3). Model I has no meaningful analytic solution if the number of individuals of any ploidy level is zero (except for the trivial solutions at *b* = 0 or *μ* = 0); in such cases, we replaced the zero relative frequency by 1/(2*N*), where *N* is the sample size for a given species (i.e. assuming that one individual of the missing ploidy level would be found if the sample size were doubled). The values obtained thus have to be viewed as upper/lower limits of true values. The model can always be solved, as in both these cases the number of parameters to be estimated equals the number of independent data points (proportion of triploids in the minimal model, and proportions of triploids and tetraploids in model I). We examined whether the values found are biologically meaningful.

In model II, parameter values were estimated using a parameter optimization procedure in which we compared observed frequencies of individual ploidy levels with the models' predictions, using as a criterion to be minimized, where *n _{i}* is the observed number of individuals of the

*i*th ploidy in a given species and

*v*is the element of the dominant eigenvector corresponding to the

_{i}*i*th ploidy of the matrix

*M*(equation (2.3)). Parameter optimization was done by a quasi-Newton method of Broyden, Fletcher, Goldfarb and Shanno (BFGS) as implemented in the function

`optim`in R v. 2.9.1 (module stats; [20]). All the parameters we work with have natural constraints: probabilities of unreduced gamete formation are bound by zero and unity; relative fertility has a lower bound of zero, and although it can theoretically have values greater than unity, in accordance with the empirical data we set its upper limit to unity as well. We therefore used a modified version of the BFGS algorithm [21], which allows box constraints by which each variable is given a lower and/or upper bound. We always examined reproducibility of the solution by starting from multiple different likely initial parameter values. We considered fit of the model to a particular dataset to be good if the optimization criterion was smaller than 0.01, poor but still informative if the optimization criterion was between 0.01–0.04, and models with even poorer fits were considered unrealistic with their resulting parameter values doubtful. Here we estimated three parameters (

*b, c*and

*μ*) in the system with four independent data points, i.e. we examined whether the models match the data (closeness of the fit) and whether the values found are biologically meaningful.

## 3. Results

### (a) Predictions of the gametic model with only diploids capable of producing unreduced gametes (model I)

The behaviour of the model shows a region of coexistence of all three cytotypes (when *μ* < 1−*b*) and a region where only tetraploids exist (figure 1). It can be shown analytically that it has a threshold at *μ* = 1−*b*, close to which equilibrium proportions of individual ploidies are extremely sensitive to the parameter values (figure 1). Equilibrium proportions of cytotypes depend on both parameters *μ* and *b*, and are primarily affected by the frequency of unreduced gametes in the diploid (figure 1). The effect of the relative fertility of tetraploids is highly nonlinear: close to the threshold, proportions of ploidies are highly sensitive to *μ*, but otherwise the effect of *μ* is rather weak (figure 1).

Triploid plants are never predicted by the model to have high frequencies. The proportion of triploids is always affected by the proportion of unreduced gametes in the diploid (*b*), and, in the variant of model I with a general value of triploid fertility (i.e. *μ*_{3} ≠ *μ*_{4}), also, but to a lesser degree, by fertility of triploids (*μ*_{3}) relative to tetraploids (*μ*_{4}). Because triploids are inherently unstable (only half of their progeny are triploid if unreduced gametes are not generated, and an even lower proportion if they are generated; see the electronic supplementary material, table S2 and figure S2), their frequency is governed primarily by the parameters of their formation, not by their own fertility.

### (b) Predictions of the gametic model with both diploids and triploids capable of producing unreduced gametes (model II)

In this scenario, there can be five ploidy levels, from diploid to hexaploid. In principle, the system has six parameters, but assuming *μ*_{3} = *μ*_{4} = *μ*_{5} = *μ*_{6} (which we will denote *μ*) leaves a system with three parameters. If relative survival rates of higher ploidy levels (*μ*) are sufficiently far from unity, the system has a stable ploidy structure, which is determined by all three parameters. In this case, either diploids or tetraploids can dominate the system; proportions of individual ploidies are largely determined by the proportion of unreduced gametes produced by diploids (*b*) and much less by the proportion of unreduced gametes from triploids (*c*; figure 2).

Similar to the previous case, there is a threshold value of the relative fertility of higher ploidies (*μ*), above which the ploidies are not able to coexist (figure 2). In contrast to the previous case, however, above the threshold the matrix *M* has two identical eigenvalues corresponding to two different stable structures: tetraploids only and hexaploids only. The incidence and proportions of these two ploidies are thus not determined by the structure of the matrix only, but also by the initial conditions. No other ploidy is able to persist above the threshold. The position of the threshold is affected only by the values of *b* (proportion of unreduced gametes in the diploid) and relative fertility of higher ploidies (*μ*), and is independent of the value of *c* (proportion of unreduced gametes in the triploid).

### (c) Fitting gametic models to the empirical data

In species with 2*x*, 3*x* and 4*x* ploidy levels, model I could be successfully fit to all species/populations in the dataset with the exception of two populations of *Galax urceolata* (table 2; see also figure 3). Aside from one other population of *Ga. urceolata*, the estimated values of *b* range from 5 × 10^{−4} to approximately 5 × 10^{−2}. Relative fertilities of higher ploidies in model I were very high, close to unity (but smaller than 1–*b*, which is the condition for coexistence of all ploidies). In species with only 2*x* and 3*x* ploidy levels, the predicted values of relative fertility of higher ploidies were considerably smaller, although typically non-zero. The minimal model (i.e. keeping the relative fertility of higher ploidies at zero), yielded higher values of *b* and equally good (if not better) fit.

When model II was fit to species/populations with 2*x*, 3*x* and 4*x* ploidies, it mostly predicted the value of *c* (proportion of unreduced gametes of the triploid) to be zero, which confirmed that model I is sufficient to explain the ploidy structure in these species. In a few cases, however, it predicted a non-zero value of *c* and consequently also a low proportion of hexaploids (which were not found in the field using our sample sizes). In *Chamerion angustifolium, Knautia arvensis* P3 and *Melampodium leucanthum,* data seem to be compatible both with model I with *c* = 0 (i.e. no unreduced gametes of the triploid) and model II with *c* > 0 (table 2). Relative fertilities of higher ploidies were predicted by model II to be close to unity in most species. Only in the two populations of *Ga. urceolata* to which model I could not be fit was the relative fertility of tetraploids estimated by model II considerably smaller than 0.9.

A satisfactory fit using model II could be obtained in nine out of 11 species/populations with 2*x*–6*x* ploidy levels (table 2). Values of *b* fell into a similar range as in the 2*x*–3*x*–4*x* species, but proportions of unreduced gametes in triploids (*c*) were (rather expectedly) predicted to be considerably higher than proportions of unreduced gametes in diploids (*b*). They ranged widely, from almost zero (in *K. arvensis* P4) to over 0.7 (in *Solidago altissima*). The fertility of tetraploids was again predicted to be rather close to unity in most of the species. *Pilosella echioides*, although its ploidy structure was predicted very well (table 2), had quite exceptional predicted parameter values (very high *b* and very low fertility of higher ploidies). *Oxycoccus palustris* and *P. echioides* P1 could not be fit by model II at all. In *Oxycoccus,* high frequency of pentaploids is not compatible with a complete absence of triploids (and hence diploids) under the simple gametic model.

## 4. Discussion

The last decades have seen several attempts to model the conditions, allowing establishment of polyploid derivates within a population of diploid/lower-ploid ancestors [8–15]. Our work adds to these efforts and presents a gametic model for different diploid-autopolyploid systems occurring under natural conditions. We show that simple assumptions of the gametic process yield a model that predicts a number of patterns observed in the field, and parameters of this process can be estimated (either analytically or using an optimization procedure) from proportions of individual ploidy levels under natural conditions. While the model involves a number of simplifications (absence of crosses between ploidy levels, and identical fertility of all higher ploidies), it does capture the key processes in autopolyploid generation and persistence, *viz*. formation of unreduced gametes and reduced fertility of higher ploidies. Its good fit to most of the data indicates that these processes are sufficient to explain ploidy heterogeneity that these species show.

### (a) Gametic model of polyploidy

Ploidy composition predicted by the gametic model is largely driven by the parameters of even ploidy levels (which are more common), and diploids in particular (which always represent the starting point). The low effect of odd ploidy parameters on overall ploidy composition as predicted by the model is caused by both the rarity of odd cytotypes and their inherent instability (i.e. the formation of gametes with different ploidy levels, and consequently, ploidy-variable progeny). This confirms earlier findings that while triploids contribute to polyploidization of a population, the dynamics of the process itself does not strongly depend on their parameters [12].

The gametic process also constrains the proportions of triploids and tetraploids given the proportion of diploids. If tetraploids are not detected and triploids present, it indicates a process that precludes tetraploids from forming or bearing seeds, such as very low survival of tetraploids, or triploid sterility [19]. If the proportion of triploids is higher, the fit necessarily decreases, as the gametic model requires some tetraploids to form. It should also be noted that in systems with no tetraploids, the value of relative fertility of higher ploidies (*μ*) depends on the upper estimate of the frequency of tetraploids, and is thus partly arbitrary. The same type of constraint applies to higher ploidy levels.

The gametic process predicts coexistence of several ploidies only if survival/fertility of higher ploidies is lower than relative production of *reduced* (i.e. normal) gametes by the diploid. Close to this threshold, the effect of the relative fertility becomes dramatic, and once crossed, higher ploidies act as a permanent sink with the whole system developing into a homogeneous population of only one higher, even polyploid level (depending on other parameter values and on initial conditions; see also [10–12]). Around this threshold, parameter values in the model are not susceptible to detection owing to the extreme sensitivity of predicted ploidy proportions to parameter values.

### (b) Fit of the model to the data

The major novelty of our study is the comparison of *in silico* results with real ploidy compositions observed under natural conditions. The validity of the model can be assessed by ascertaining whether (i) predicted ploidy frequencies match with the empirical ones, and (ii) the inferred parameter values are compatible with values obtained using more direct (experimental) approaches.

With a few exceptions (*Ga. urceolata*, *O. palustris* and *P. echioides*, which often showed low fit), the predicted frequencies of unreduced gamete production for diploid plants (5 × 10^{−4} − 5 × 10^{−2}) matched the empirical data well. The mean production of diploid gametes across non-hybrid flowering plants is 5.6 × 10^{−3} [6], which is not very far from the values predicted by our models (mean 1.66 × 10^{−2}, median 8.9 × 10^{−3}, across all plants/populations and models with good fit; the unrealistically high values for the above-mentioned species excluded). Under natural conditions, the observed frequencies of diploid gametes vary substantially among species as well as among conspecific individuals. For example, the proportion of unreduced pollen in individual clones of *Vaccinium* sect. *Cyanococcus* varied between 0.01 and 0.28 [22], while mean population-level frequencies of unreduced pollen in *Achillea borealis* ranged from 3 × 10^{−4} to 5.4 × 10^{−3} but were as high as 0.158 in individual plants [23].

Conversely, external knowledge of the reproductive and fitness parameters (such as the proportions of unreduced gametes and survival rates of higher ploidy plants) can be used as a further test of the model realism. In such a case, predicted frequencies of minority cytotypes follow directly from substituting these parameters into the model and can be compared with the field-observed values.

Sample size permitting, the model was fit both to species (i.e. cumulative ploidy frequencies across all populations) and individual populations. Despite the greater effect of stochastic processes (which can cause the elimination of a particular cytotype and preclude model fitting; 3*x* and 4*x* populations of *Ga. urceolata* P2 may be an example), the latter show slightly better fit perhaps owing to more accurate field data on ploidy composition at finer spatial scales (cf. population and species data for *Gymnadenia conopsea* or *K. arvensis*). In a few cases, some species/populations could not be fit with biologically meaningful parameters, indicating incompatibility between what is occurring in nature and one or more model assumptions. First, individual cytotypes, especially in mixed-ploidy populations and/or in contact zones, are not necessarily fully reproductively isolated and can hybridize. Whereas there often is a strong breeding barrier between diploids and tetraploids (i.e. the so-called triploid block; [19]), increase in genome copy number can be associated with relaxation of intercytotype reproductive barriers [24]. Consequently, interploidy interactions can play more important roles in shaping ploidy composition than do gamete-driven processes. *Oxycoccus palustris* [25] and *P. echioides* [26] probably fall into this category and both show poor fit or unrealistic predicted parameter values. Moreover, the model predicts, in agreement with most field data, large prevalence of even ploidy levels. Therefore, fit was usually unsatisfactory (or yielded unrealistic predicted parameter values) in groups with a high proportion of odd ploidy levels (*Ga. urceolata*, *P. echioides*). Additionally, the current formulation always predicts a non-zero proportion of diploids in species with ploidy heterogeneity, therefore entailing poor or even no fit in species/populations lacking diploids (e.g. *O. palustris*). In general, poor fit to the gametic model can be used as a strong indication of action of other processes (hybridization, strongly unequal fertility or survival of higher ploidy levels, etc.), in the species in question.

### (c) Implications

The main feature of the gametic model is the requirement that higher ploidies must be less fertile relative to diploids if all ploidies are to be maintained in the system (see [12]). While this represents a difficulty for model fitting, it is likely to be a real process occurring in the field and shaping the evolution of polyploids. The relative fertility of polyploids needs not be very low (it only must be lower than the proportion of *reduced* gametes in the diploid), and if this condition is met, its effect on the proportions of all ploidies is relatively mild. The threshold effect of relative fertility of higher ploidies seems to be a universal property of all models, if immediate reproductive isolation of each ploidy level is assumed.

After an (auto)polyploidization event, it is likely that the newly formed polyploids will be less fertile than diploids [27]. However, individual selection on tetraploids will act to increase this fertility in polyploids, driving the value to the threshold; once the threshold is reached, the whole population will suddenly switch to a high even polyploid level, leaving no diploids. While we have no field data to support the existence of this phenomenon of a sudden transition, it is a universal feature of gametic models (see [10]). We believe that it may indeed account for the steady increase of ploidy levels during the evolution of flowering plants (see [14]).

Our model does not support the widely held assumption that higher odd ploidy levels ought to be less fertile than their even ploidy counterparts. Our results show that the saw-toothed frequency distribution of ploidies can be generated by simple ploidy-specific differences in the frequency of gametic types in otherwise equally fertile odd and even polyploids. Observed predominance of even ploidy plants observed under natural conditions can thus simply be owing to inherent instability of the gametic process in their odd ploidy counterparts (although in several well-documented cases the triploid fertility is indeed much lower [7]).

Although the model relies on some important assumptions (e.g. complete reproductive isolation among different cytotypes), it seems likely that the conditions under which autopolyploids can become established and coexist are less restrictive than previously assumed [28]. This can at least partly explain the prevalence of polyploids and fairly frequent occurrence of mixed-ploidy populations in natural conditions. The lower fit of the model to *in situ* ploidy composition in some plant groups can serve as a good indication of more complex ploidy-generating processes (including interploidy hybridization) involved. As the role of such processes in the evolutionary dynamics of mixed-ploidy species/populations needs to be better understood, a simple gametic model may also be used to identify such cases from field-collected ploidy frequency data.

## Acknowledgements

The support was provided by the Academy of Science of the Czech Republic (long-term research development project no. RVO 67985939), institutional resources of Ministry of Education, Youth and Sports of the Czech Republic for the support of science and research, and the Czech Science Foundation (project 206/09/0843). The William J. Fulbright Commission and University of Michigan supported T. H. during his sabbatical stay at the Department of Ecology and Evolutionary Biology. We thank all our colleagues for permission to use their unpublished data on ploidy frequencies (Martin Hanzl, *Knautia*; Jana Krejčíková, *Campanula*, *Oxalis*; Karol Marhold, *Cardamine*; Petr Vít, *Cerastium*). Jonathan Rosenthal kindly improved the English. We thank two anonymous referees for their comments that helped us to clarify the paper.

- Received October 9, 2012.
- Accepted November 2, 2012.

- © 2012 The Author(s) Published by the Royal Society. All rights reserved.