Clones on the run: The genomics of a recently expanded partially clonal species

Why species that in their core areas mainly reproduce sexually become enriched with clones in marginal populations (“geographic parthenogenesis”) remains unclear. Earlier hypotheses have emphasized that selection might promote clonality because it protects locally adapted genotypes. On the other hand, it also hampers recombination and adaptation to changing conditions. The aim of the present study was to investigate the early stages of range expansion in a partially clonal species and what drives an increase in cloning during such expansion. We used genome‐wide sequencing to investigate the origin and evolution of large clones formed in a macroalgal species (Fucus vesiculosus) during a recent expansion into the postglacial Baltic Sea. We found low but persistent clonality in core populations, while at range margins, large dominant clonal lineages had evolved repeatedly from different sexual populations. A range expansion model showed that even when asexual recruitment is less favourable than sexual recruitment in core populations, repeated bottlenecks at the expansion front can establish a genetically eroded clonal wave that spreads ahead of a sexual wave into the new area. Genetic variation decreases by drift following repeated bottlenecks at the expansion front. This results in the emerging clones having low expected heterozygosity, which corroborated our empirical observations. We conclude that Baker's Law (clones being favoured by uniparental reproductive assurance in new areas) can play an important role during range expansion in partially clonal species, resulting in a complex spatiotemporal mosaic of clonal and sexual lineages that might persist during thousands of generations.

Alternatively, clonal dominance can result from the advantage of reproducing without mates during establishment in previously noncolonized areas (Baker's Law) (Baker, 1955;Pannell et al., 2015;Tilquin & Kokko, 2016). Such clone dominance has been previously found through range expansion modelling, where partially clonal species with low rates of migration form a clonal wave that expands ahead of a sexual wave, establishing large patches of clones that prevail for thousands of generations at range margins (Rafajlović et al., 2017). Importantly, stochastic processes at the expansion front (Excoffier et al., 2009;Slatkin & Excoffier, 2012), the balance between recombination and clonality (Eriksson & Rafajlović, 2021;Peck et al., 1998;Vrijenhoek & Parker, 2009), and different dispersal strategies (Ibrahim et al., 1996;Paulose & Hallatschek, 2020;Rafajlović et al., 2017) will influence the population genetics of such expansions of partially clonal species. Mathematical modelling can be used to illustrate the interplay between these evolutionary and demographic processes and support the interpretation of empirical population genomic data. Likewise, the inbreeding index F IS , observed (H o ) and expected (H e ) heterozygosity, and the linkage disequilibrium are straight-forward indices that estimate the levels of clonality and genetic variation. In a large, well-mixed sexual population, H e and H o are expected to be approximately equal and F IS will be close to zero, while negative values of F IS will indicate the presence of cloning (Arnaud-Haond et al., 2020;Balloux et al., 2003;Stoeckel & Masson, 2014). Thus, analysing the geographic patterns of F IS , H o and H e in areas of recent range expansions, in combination with range expansion models, will allow us to elucidate the evolutionary history of a colonization event, specifically tracing the origin and spread of clones contributing to patterns of geographic parthenogenesis.
The postglacial Baltic Sea is one example of a marginal area where several common species of seaweed and macroalgae have switched from a dominance of sexual reproduction outside the Baltic Sea to the occurrence of large clones inside (Johannesson & André, 2006). One outstanding example is the brown algae Fucus vesiculosus, a common intertidal macroalga that expanded from the Atlantic into the postglacial Baltic Sea less than 8000 years ago and switched from seemingly fully sexual outside to partially clonal reproduction inside the Baltic Sea (Pereyra et al., 2009;Tatarenkov et al., 2005). This species is dioecious (separate female and male individuals) with external fertilization. The zygotes formed have restricted dispersal (a few meters, Serrão et al., 1997), but occasionally, whole individuals drift long distances (Rothäusler et al., 2015).
The distribution is essentially one-dimensional, in the Atlantic being restricted to the lower intertidal and in the Baltic Sea to the shallow subtidal. However, populations are usually dense and continuous along rocky shores. An earlier study of populations from the southern Baltic Sea and outside the entrance did not find evidence of inbreeding (Tatarenkov et al., 2007). Furthermore, earlier studies using microsatellite markers have identified a complex pattern of sexually dominant, clonal dominant or mixed populations but failed to explain the evolutionary causes and consequences of this obvious case of geographic parthenogenesis (Ardehed et al., 2016;Pereyra et al., 2013). Here, we combine extensive SNP-genotyping of a large number of Baltic and nearby Atlantic populations to address these issues with more powerful tools. To support the interpretation of our empirical data, we also expanded a mathematical model (Rafajlović et al., 2017) earlier used to examine range expansion of partially clonal species. Both the new and the old models predicted our empirical observations in F. vesiculosus, and strongly supported the role of Baker's Law and stochastic processes in shaping the range expansion and forming compelling patterns of geographic parthenogenesis in this partially clonal species.

| DNA extraction and genotyping
Using 2b-RAD sequencing (Wang et al., 2012), we successfully genotyped 962 individuals (896 plus technical replicates) from 47 sites representing the Baltic Sea distribution and nearby areas of the Atlantic Ocean (Roscoff, Bangor, southern Norway and further into the most northern and eastern distributions in the Baltic Sea) (Table S1). Genomic DNA was extracted from 20 individuals per locality (Table S1) and re-extracted from four of those 20 individuals to use as technical replicates (i.e. replicated extraction, library preparation and sequencing). We followed Panova et al. (2016) for genomic DNA extractions. Thus, we ground dried algal tissue and washed it twice with 100% acetone before digestion in CTAB buffer. This was followed by extractions using the DNA Plant II Extraction Kit (Qiagen) with an extra cleaning step using the DNA Clean and Concentrator kit −25 (Zymo). We assessed the DNA integrity by electrophoresis in a 1% agarose gel and the quality using Nanodrop.
We used the computer cluster 'Albiorix' at the University of Gothenburg, Sweden, for all our bioinformatic analysis, including a modified pipeline for trimming, filtering and genotyping, available at https://github.com/crust acean a/TheFu cusPr oject and originally written by Mikhail Matz (https://github.com/z0on/2bRAD_denov o/blob/maste r/2bRAD_README.sh). Sequences were trimmed and quality filtered before being mapped to an unpublished F. vesiculosus draft genome assembly (NCBI Assembly Number ASM1454947v1). This genome has been used to map F. vesiculosus RAD-seq reads with high success (Kinnby et al., 2020) despite the assembly's relatively low completeness and fragmentation (assembly statistics available at NCBI). Considering the potential impact on the SNPs yield, we performed a preliminary de novo 2bRAD analysis for a quick comparison of the number of SNPs with the reference-based one, but the de novo approach yielded less SNPs than the reference-based analysis, so the latter was kept for subsequent steps. We subsequently called the genotypes with the GATK v. 4.1.0.0 pipeline and used the technical replicates to extract 100% reproducible SNPs across all replicates, with heterozygotes in at least two replicate pairs. We used 18,866 SNPs as 'true' SNP dataset to train the error model for non-parametric quantile-based recalibration of variants.
The recalibration rendered 325,885 SNPs, from which loci with poor coverage were removed, while retaining only bi-allelic polymorphic sites and loci genotyped in ≥90% of all individuals. Further filtering excluded loci and individuals with >5% of missing data using the R package "radiator" v. 1.1.2 (Gosselin et al., 2020) and rendered a final set of 15,409 SNPs. The genotyping accuracy based on the overall match between replicates was 96.56%, that we used as a preliminary threshold to identify clones. After variant recalibration, we removed the technical replicates and thinned the dataset to retain one SNP per RAD fragment. All raw 2b-RAD sequences are deposited in the Sequence Read Archive Repository at the National Centre for Biotechnology and Information (NCBI BioProject Number PRJNA629489).

| Evolution of sexual and clonal lineages
We used a genetic distance-based method to infer the genetic distance threshold separating sexually and clonally recruited genotypes (Kamvar et al., 2014). We first determined the number of clones in our dataset and the geographic patterns of clonal expansion with a complete dataset comprising 962 individuals distributed across all 47 localities. We initially calculated genetic distances between individuals using the function diss.dist in the R package "Poppr" v.2.9.3 (Kamvar et al., 2015), which employs Prevosti's genetic distance to estimate the fraction of different SNPs between samples.
We used this genetic distance matrix to plot a Neighbour-joining (NJ) tree to identify the geographic distribution of clonal lineages and clonal clusters. We then calculated a conservative maximum genetic distance as cut-off for assigning clones using the 'cutoff_predictor' function in Poppr. The 'cutoff_predictor' function will filter out genotypes by degree of similarity using three clustering algorithms (Nearest Neighbour, UPGMA and Farthest Neighbour) to search the top fraction of genetic differences and find the largest genetic difference between MLGs. Our calculations showed consistent values using Nearest Neighbour and UPGMA clustering, and we used these to calculate the average between the distances spanning that largest genetic distance and to estimate the cut-off threshold defining the clonal lineage threshold (Kamvar et al., 2015). By including technical replicates in this estimation, it allowed us to account for sequencing errors. After estimating the cut-off threshold, we identified similar multi-locus genotypes (MLGs) using the mlg.filter function from Poppr to plot the frequency distribution of distances between individuals in relation to the estimated threshold (0.0197; supplementary Figure S1). (This threshold resulted in 792 contracted MLGs in 896 individuals, after removal of technical replicates.) The rationale for using this method to detect clonality is that if the distribution of genetic distances among individuals shows distinct peaks towards low distances instead of a strict unimodal distribution, these peaks are likely to reveal somatic mutations resulting in low distances among slightly distinct MLGs deriving from the same reproductive event (Arnaud-Haond et al., 2007;Douhovnikoff & Dodd, 2003;Meirmans & Van Tienderen, 2004).
We further corroborated the presence of clonality by testing each population for linkage disequilibrium using the standardized index of association (r d ) (Agapow & Burt, 2001;Brown et al., 1980). Doing this, we defined all individuals from the same locality as one population, but treated the six clonal clusters as separate populations. Random sampling of the genome is expected to have a mean index of association close to zero for panmictic populations as opposed to populations under clonal reproduction (Kamvar et al., 2015).
However, genome-wide sampling of potentially thousands of anonymous markers will likely create deviations from zero (random mating) even within sexual populations. Nonetheless, genome-wide linkage disequilibrium is expected to be considerably lower in sexual populations relative to clonal lineages.
The NJ tree revealed separate clonal clusters within the "Baltic" clade suggesting the repeated independent evolution of clonal lineages. However, the tree could not reveal reticulations among these clusters and therefore we plotted a Minimum Spanning Network (MSN) with "Poppr" using our genetic distance matrix, to reveal inter-and intraclonal reticulations and whether or not the clonal lineages had evolved independently. We additionally characterized the genetic variation and clonality by calculating the relative levels of observed (H o ) and expected (H e ) heterozygosity, and inbreeding coefficient (F IS ) for each population using the R package "radiator".
(Please note that the estimates of H o and H e are not showing absolute levels, because the RAD dataset had been filtered for polymorphic sites and with a minimum threshold of minor allele frequency).

| Statistical analysis of population genetic estimators
We performed statistical tests for overall differences in genetic variation (H e ) and observed heterozygosity (H o ) among the three geographic areas, NE Atlantic, Transition zone, and Baltic Sea both including and excluding the clonal lineages using the GenoDive v.3.05 software (Meirmans, 2020). This included permutation tests to identify differences among geographic groups. OSx-statistic (Goudet, 1995) was used as a test statistic comparing the sum of the squared differences over all pairwise combinations of groups.
Permutations took place by randomizing the populations over the groups.

| Spatial and temporal modelling of clonal expansion
We explored the dynamic relationships between H e , H o and F IS for a modelled partially clonal species expanding over a one-dimensional habitat (a coastline) of 500 patches. Our model was a development of the model used in Rafajlović et al. (2017) and it differed in the following aspects. First, we increased patch size from N = 1 to N = 100 individuals maximum resulting in a larger meta-population consisting of well-mixed local populations. Second, we included an old 'core population', the first 50 patches, from which the expansion started instead of starting from 100 individuals with alternating sexes.
Third, we used a linear, instead of a circular row of patches, because a linear expansion in one direction was appropriate to reflect the early phases of range expansion. Fourth, we added a possibility of a successive loss of sexual function in clonal lineages. Fifth, in the earlier model we only followed the sex ratio, but now we in addition followed the temporal and spatial genetic patterns of expected heterozygosity (H e ), observed heterozygosity (H o ), and the inbreeding index F IS by explicitly modelling the evolution of allele frequencies at L = 100 neutral, bi-allelic and freely recombining loci. Finally, we also introduced mutations at each locus, both in clonally and sexually produced offspring following a symmetric two-allele model, with a per allele, per offspring, per reproductive season probability fixed to a small value (μ = 0.000025) in all runs.
The model parameterisation was guided by life-history characteristics of F. vesiculosus, including separate sexes, larger perindividual investment in sexual, rather than in asexual, reproduction (unless the sexual function was partially or completely lost), dispersal mostly to nearby patches but with occasional long-distance dispersal, and a perennial lifespan (Table 1).

| Model starting conditions
We started all simulations with a burn-in period that we initiated with setting the expected heterozygosity at each locus to H e,0 = 0.2 , and applying Hardy-Weinberg and linkage equilibria. To generate the starting genotypes at each locus, we chose uniformly at random among the two alleles (0 and 1) with probability long-range dispersal occurred with probability p L , and short-range dispersal with probability 1 − p L . Short-range dispersal was always to a neighbouring patch, and both directions were equally likely. Longrange dispersal, equally likely in either direction, was at a distance randomized from an exponential distribution with mean assumed to be much larger than the distance between neighbouring patches, TA B L E 1 Model parameters with notations used, brief explanations and values explored.
Starting per-locus expected heterozygosity in the source patches, i = 1, 2, -K 0 0.2 but much smaller than the distance across the complete habitat, and we set δ = 50 throughout. Because we modelled a discrete habitat space, we discretised the exponential distribution of long-range dispersal distances similarly as in Eriksson and Rafajlović (2021) and long-range dispersal was required to be at least double the distance of short-range dispersal and less than the size of the habitat (this was achieved by truncating the probability distribution at distance two, and K from below, and above, respectively). Individuals that were migrating away from the habitat boundaries were placed in the nearest boundary patch. Furthermore, immigrants could only inhabit a patch they arrived at if there was space available, which was regulated by random removal of excess immigrants.

| Sexual and clonal reproduction
denote the sexual function of male k, the number of males and of females, respectively, in patch i in reproductive season t. Here, index i = 1, 2, … , K is the patch number. When the sexual function of all males was maximal, that is, when , a randomly chosen female experienced, on average, f (t) encounters, and this number was bounded from above by S. Note that, due to the stochastic nature of encounters (governed by males' sexual function and random dispersal of sperm) and of fertilization (governed by females' sexual function), it could happen, by chance, that a female with a smaller sexual function gave rise to more sexually produced offspring than a female with a larger sexual function. Conversely, the probability that one encounter resulted in fertilization was where s (i,j) f,s (t) is the sexual function of female j in patch i in reproductive season t. When the sexual function of all females was maximal, that is, when , the probability that an encounter resulted in fertilization was To trace the number of clonal lineages in a patch, we assigned an ID number to each individual, with each sexually produced offspring receiving a new ID. We assumed that sexually produced offspring were either males or females with equal probability, and that the sexual function of all sexually produced male and female offspring was maximal s m,s = s f,s = 1. (This is a technical simplification of the model but it will not promote clonal propagation.) In addition, individuals were also capable of reproducing clonally, and in a given reproductive season, each individual also contributed with C clonally produced offspring which inherited the sex, genotype and ID of their parent. However, a clonally produced offspring could differ from the parental genotype due to a new mutation (see We were interested in parameter combinations under which clonal lineages form, and one requirement for this was low rates of dispersal (Rafajlović et al., 2017), and hence, we present model results for a subset of parameter values that resulted in the formation and establishment of one or more local clonal lineage that persisted during at least 500 reproductive seasons.

| The genomic landscape of Fucus vesiculosus in a marginal environment
Overall Five clusters were restricted to a single locality, while cluster 45 and its largest clonal lineage, the most distal in the NJ tree, was spread over six localities (45a-45f in Figure 1b, and see Figure S6) and >600 km of coastline. Despite restricted sampling in each site, single pairs of MLGs belonging to the same clonal lineage were also found in the southern Baltic Sea, the Transition zone, and in one of the Atlantic sites ( Figure 1a) indicating that clonality also occurs in core areas but in low proportion.

| Spatial patterns of H e , H o F IS and linkage disequilibrium
In sexual lineages, both H e and H o tended to decrease along the colonization route from the Atlantic to the inner parts of the Baltic Sea, but the differences between Atlantic and Baltic sites were minor while the Transition zone sites peaked somewhat (Figure 2,   (Table S3). Clonal cluster 40 showed considerably lower r d than the other clonal clusters but still significant.

| Model outcome with short-range dispersal
Under restricted and short-range dispersal, the simulated expansion of a partially clonal species over an empty one-dimensional habitat showed two main results: (i) large clonal lineages remained absent in the core population (patches 1-50) under the model assumption of a demographic advantage of sexual reproduction, and (ii) a clonal wave formed and spread ahead of a sexual wave and became dominant during, and long after, 10,000 reproductive seasons (Figure 3ac). Eventually, after >>10,000 reproductive seasons, we expect that the subsequent spread of sex will assimilate the clonal wave, as partners for sexual reproduction slowly become available. Note that, during the simulated 10,000 reproductive seasons, the blue area (sexual wave with fully or partly sexual reproduction) increases slowly in size in Figure 3a

| Adding occasional long-range dispersal
With occasional long-range dispersal, the initial clonal wave spread to new patches and occasionally additional clonal waves originated from encounters of sexual mates introduced by long-range dispersal (Figure 4 and Figure S8)  Figure 1b. If the same site includes one clonal cluster and one sexual population, or two clonal clusters, these have different numbers, and note that sample 45 was present in six of the sites (see Table S1). Green points represent dominant clonal clusters, red dashed lines indicate overall averages of H e and H o (clonal clusters not included) and F IS = 0, respectively. Error bars in H o illustrate SD and in F IS represent variance among the 15,409 2bRAD loci in each population (where clonal clusters are treated separately when sampled from the same sites as sexual individuals).  Figure S8 for H e , F IS , and the evolution of sexual function). Captions are the same as in Figure 3. The rows differ by the total migration rate (m T ) and the probability that a migrant disperses by long-range dispersal (p L ) as depicted in the figure. The columns differ by the rate of loss of sex (s Loss ) as indicated in the figure. The mean long-range dispersal distance was set to 50 patches. The remaining parameters are as in Figure 3.
rate would be very long and the clonal lineages would therefore appear permanent (Figure 4c). Under occasional long-range dispersal, heterozygosity in the core population was also sensitive to rates of short-distance dispersal (Figure 4g-l) and under most scenarios variation was already lost upon the formation of clonal waves. However, stochasticity played an important role, not least for higher rates of short-distance dispersal, and occasionally clonal lineages were formed with relatively high H o (Figure 4j-l). As in the other simulations, over time, persistent clonal lineages accumulated new variation from mutations. Under occasional long-range dispersal, patterns of individual runs showed similar trends as with only short-range dispersal, with an initial strong drop in genetic variation before one clonal lineage was formed followed by additional ones later, although the overall structures were highly complex (Figures S11-S13).

| DISCUSS ION
Baker's Law states that selfing and partially clonal species are more successful colonizers than species with entirely sexual reproductive strategies due to their uniparental advantage (Baker, 1955;Pannell et al., 2015;Tilquin & Kokko, 2016). Modelling has shown that the mechanism of Baker's Law is also relevant during the range expansion of partially clonal species (Rafajlović et al., 2017). Despite cloning being less successful than sexual reproduction in core areas, a single clonal lineage will form during the early phase of the expansion if dispersal is restricted and only stepwise. This clonal lineage will fill up most of the new territory forming a clonal wave before it is later caught up and assimilated by the spread of mates and the formation of a sexual wave. Adding occasional long-range dispersal provides opportunities for formation of multiple clonal lineages.
Moreover, a clonal lineage that loses its sexual function before being caught up by the sexual wave will resist assimilation and persist, until being outcompeted by lineages that are demographically or evolutionary more advantageous.
Does the principle of Baker's Law better explain the pattern of geographic parthenogenesis found in Fucus vesiculosus than alternative models? For example, cloning might be more favourable than sexual recruitment in extreme environments by reducing swamping of non-optimal alleles from core areas (Peck et al., 1998). However, this fails to explain our empirical observations of sexually recruited Fucus populations living side-by-side with clonal lineages at range margins. Alternatively, populations of multiple, highly heterozygous clonal lineages might be more successful than sexual populations in heterogeneous environments (Vrijenhoek, 1979). Nevertheless, observed heterozygosity in Baltic Fucus is similar in clonal and sexual lineages. On the other hand, the wide distribution of clonal cluster "45" suggests a General Purpose Genotype (Vrijenhoek & Parker, 2009), although a generalist genotype would be most successful under high dispersal (Tilquin & Kokko, 2016) while low dispersal and locally adapted populations typically characterize F. vesiculosus (Kinnby et al., 2020;Serrão et al., 1997).
When parameterised with demographic traits of F. vesiculosusthat is, predominant short-range dispersal (Serrão et al., 1997;Tatarenkov et al., 2007) occasional long-range transport (Rothäusler et al., 2015) and a general capacity of clonal reproduction-our  , only <3% putatively sexually recruited individuals were found in a natural site where this clonal lineage lives together with a male clonal lineage (Ardehed et al., 2015).
Further observations of sites where two clonal lineages of opposite sex, or clonal and sexual lineages were mixed (e.g. sites 35, 39, 41 in Figure 1b) suggest that the widespread clonal lineages somehow evade assimilation. Earlier studies suggest that the low salinity (<6‰) in the marginal Baltic areas has negative impact on sexual function, and increased rates of polyspermy have been found in sites at the Swedish coast (Serrão et al., 1999). However, sexual populations are present down to salinities of ca. 2-3‰ in Gulf of Finland (Ardehed et al., 2016), and thus low salinity does not completely inhibit sex, although it might prolong the persistence of clonal lineages. Finally, the demographic advantage of the sexually functional individuals, explaining their dominance outside the Baltic Sea, is likely to eventually outcompete the clonal lineages, albeit at a much slower pace than the assimilation process involving recombination.
Once formed, a clonal lineage essentially freezes the level of observed heterozygosity. Furthermore, in a partially clonal population the approach towards Hardy-Weinberg equilibrium is strongly reduced, and non-zero values of F IS will appear also in populations with moderate levels of clonality (Reichel et al., 2016). In addition, genetic variation (H e ) might be lost when a new clonal lineage is formed, but can later be restored by genotypic drift (Bengtsson, 2003;Reichel et al., 2016;Stoeckel et al., 2021). Our modelling results also suggest that during range expansion, clonal lineages originate from sexual populations that have already lost much of their original genetic variation due to repeated bottlenecks at the expansion front. In our empirical data, we observed a considerable loss of genetic variation (H e ) in four out of six clonal clusters. We believe that this loss is most likely explained, as in our model, by repeated bottlenecks at the front of the expansion prior to, and during, the formation of the clonal waves. Contemporary Baltic populations of F. vesiculosus are dense and F I G U R E 5 Population genetic structure with only short dispersal after 2000 reproductive seasons. Expected heterozygosity (a-f), observed heterozygosity (g-l), and F IS (m-r) were obtained 2000 reproductive seasons after the start of expansion in the model where dispersal occurs only to the nearest neighbours. The black symbols show the values averaged over 200 independent runs, the blue and red bars span between the 5th and 95th percentiles of the realized distributions in the sexual and clonal areas, respectively. Here, a clonal area was any patch in which genotypic identity based on individuals' IDs was >0.95 in at least 50% of independent model runs, conditional on that the patch was occupied in the given reproductive season. For comparison, the green symbols and error bars show the corresponding results at the end of the burn-in simulations wherein only patches 1,2,…, 50 were fully occupied and the remainder of the habitat was empty. The rows differ by the value of the migration rate, and the columns differ by the rate of loss of sex, as depicted in the figure. The remaining parameters are as in Figure 3.
continuously distributed on shallow subtidal rocky substrates, and the current genetic variation of Baltic sexual lineages likely reflects the drift-migration balance of much larger populations than during the expansion, thousands of generations ago. In contrast, the low genetic variation of the clonal lineages unaffected by gene flow reflects the original level of variation at the expansion front.
Earlier authors have highlighted that genetic processes during postglacial colonization and range expansion set the stage for contemporary genetic structures of sexual species (Excoffier et al., 2009;Hewitt, 1999Hewitt, , 2011. Here, we additionally find that in partially clonal species, range expansion can lead to the formation of clonal and sexual waves that generate a mosaic genetic landscape that may remain for thousands of generations, the details of which are difficult to predict at specific points in space or time because of its stochastic and transient nature. We also conclude that patterns of geographic parthenogenesis can be established by stochastic processes and the mechanism of Baker's Law.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors have no conflicts of interest to declare.