Ecological success of sexual and asexual reproductive strategies invading an environmentally unstable habitat

Many aspects of sexual and asexual reproduction have been studied empirically and theoretically. The differences between sexual and asexual reproduction within a species often lead to a biased geographical distribution of individuals with different reproductive strategies. While sexuals are more abundant in the core habitat, asexuals are often found in marginal habitats along the edge of the species distribution. This pattern, called geographic parthenogenesis, has been observed in many species but the mechanisms reponsible for generating it are poorly known. We used a quantitative approach using a metapopulation model to explore the ecological processes that can lead to geographic parthenogenesis and the invasion of new habitats by different reproductive strategies. We analyzed the Allee effect on sexual populations and the population sensitivity to environmental stress during the invasion of a marginal, unstable habitat to demonstrate that a complex interaction between the Allee effect, sensitivity to environmental stress and the environmental conditions can determine the relative success of competing reproductive strategies during the initial invasion and longterm establishment in the marginal habitat. We discuss our results in the light of previous empirical and theoretical studies. Author Summary Individuals can reproduce with or without sex. Very often, closely related species are distributed in a such a way that the sexually reproducing species is most frequently found in the core habitat while the asexually reproducing species is found on the edge of the habitat range. This biased distribution of reproductive strategies across a habitat range is called geographic parthenogenesis and has been observed in several species. While many processes have been proposed to explain such a pattern, a quantitative approach of the ecological processes was absent. We investigated important differences between sexual and asexual reproduction and how these differences affect the success of sexuals and asexuals invading a marginal, unstable environment. We showed that the relative frequency of each reproductive strategy in the marginal habitat depends on how much sexuals rely on population density to reproduce and how much asexuals are affected by environmental stress relative to sexuals. Our study presents a quantitative ecological explanation for geographic parthenogenesis and provides the conditions under which different distribution patterns can emerge.


Introduction
The magnitude of the Allee effect, which can be particularly important in sexual populations 117 [44; 45; 46], may play an important role in leading to geographic parthenogenesis but has been 118 essentially ignored. And finally, population sensitivity to environmental variation in marginal 119 habitats may determine whether habitat range expansion is possible and which reproductive 120 strategies can be the most successful during range expansion, and yet these factors have been 121 overlooked in particular from the empirical perspective. Because of the necessity to address 122 these processes as leading causes of geographic parthenogenesis, we developed a metapopula-123 tion model of competition between different reproductive strategies and assessed the ecological 124 success of each competing strategy during short-term invasions and long-term establishment 125 in unstable habitats, under different environmental conditions. Our model considers both the 126 Allee effect (weak and strong) that affects sexually reproducing populations and the population 127 sensitivity to stressful environmental conditions, which can be affected by genetic diversity. We 128 use the temporal mean population size during initial invasion and long-term establishment as 129 a measure of the ecological success of each competing strategy. 130 We model population dynamics in a metapopulation consisting of two qualitatively different 132 patches. One patch (South patch) is assumed to host the larger, ancestral population, while the 133 other patch (North patch or marginal habitat) is assumed to be a habitat that has been recently 134 made available to the species under focus. The availability of a patch to a certain species can 135 be affected by several factors (e.g., global warming may make northern ecosystems available 136 to species that are typical from temperate climates). In our model, the South and North not need to find a mating partner in order to reproduce, but unlike strategy A, parents 150 and offspring are genetically different (to some extent) due to recombination during ga-151 mete production (limited but existent genetic variation). Because of the higher genetic 152 variation, strategy B is assumed to be less sensitive to changes in environmental conditions 153 than strategy A. 154 (C) Obligate sexual reproduction: Individuals need to find a compatible partner to mate with 155 and parents and offspring are genetically different. Because of the difficulty in finding 156 compatible mating partners when the population density is low, strategy C is affected by Figure 1: Graphical overview of the model design. Population dynamics include strategy-specific growth, transition and migration between an ancestral population (South patch) and a marginal habitat (North patch) where environmental conditions are unstable relative to the South patch. The panels to the right represent examples of three environmental regimes to which the North patch is exposed. Five different reproductive strategies were considered (bottom rectangle): obligate apomictic parthenogenesis or clonal reproduction (strategy A), obligate automictic parthenogenesis (strategy B), obligate sexual reproduction (strategy C), facultative apomictic parthenogenesis (strategy D) and facultative automictic parthenogenesis (strategy E). Complete description of the reproductive strategies in the main text. The subscripts/superscripts S and N indicate variables characterizing the South and North patches, respectively. The South patch differs from the North patch in its carrying capacity (K S > K N ), probability of occurrence of stressful events (p S = 0 and p N > 0), maximum environmental stress level (v S = 0 and v N > 0) and effective stress level (∆ S = 0 and ∆ S ≥ 0). M represents migration between patches.
(D) Facultative apomictic parthenogenesis: Individuals can transition between strategies A  It might be argued that strategy B (automictic parthenogenesis) is more sensitive to environ-170 mental stress than strategy A (apomictic parthenogenesis) because of the loss of heterozygosity 171 in strategy B during reproduction. However, we consider the case where genomic diversity 172 produced by automictic parthenogenetic reproduction is more beneficial than genomic homogeneity under natural selection caused by environmental stress because homozygosity caused 174 by automictic parthenogenesis introduces new phenotypes into the population.

175
Different strategies can affect one another indirectly through their effects on the total popu-176 lation density and transitions between reproductive strategies (sexual to asexual and vice-versa), 177 as explained below. For ease of reference, we define subpopulation as the fraction of the patch-178 specific total population that is composed of individuals with a specific reproductive strategy.

179
Population dynamics 180 Population dynamics follow a logistic population growth model (Ricker model) with a variable 181 growth rate, which is affected by environmental quality (carrying capacity), the Allee effect on 182 strategies C-E, and environmental effects. Additionally, the total population density is affected 183 by migration from/to the South patch. Change in subpopulation size from generation t to 184 generation t + 1 due to growth is defined by the following difference equation: where the exponential term determines the logistic population growth rate based on the Ricker 186 model and is equivalent to λ in exponential growth models, M s p i ,t is the net change in population 187 size due to migration South-North or vice-versa and T s p i ,t is the number of transitions from 188 sexual to asexual reproduction and vice-versa in strategies D-E. In strategies A-C, T s p i ,t = 0. such that each individual can produce on average at most λ = e 1.1 ≈ 3 offspring per generation 208 ( Figure S1).  The Allee-dependent growth rate (A s p i ,t ) is lowest (r min ) when the number of sexuals N + p i ,t → 216 0 and increases as N + p i ,t increases, until it reaches a biological limit (r max ). This dynamic growth 217 rate is defined by the following equation: where α is the curvature of the function and β · K p i is the population size (relative to the 219 carrying capacity K p i ) at which A s p i ,t = r max , that is, when growth rate reaches its biological

226
Our Allee model of growth rate has the following properties: (i) r s p i ,t = r min when N + p i ,t = 0;

227
(ii) r s p i ,t = r max when N + p i ,t = β · K p i ; and (iii) r s p i ,t = 0 (no growth) when We set r min = −0.1 (such that λ = e −0.1 ≈ 0.9) and assume that a very small subpopulation where p p i is the patch-specific probability of deviation from optimal environmental condi- low in strategy C (sexual reproduction; φ C = 0.5) as a consequence of their genetic/phenotypic 269 diversity. Note that strategy B differs from strategy A only in its sensitivity (φ B < φ A ). This 270 sensitivity effect is assumed to be a linear function of the environmental stress level, such that: In facultative parthenogenetic strategies (D and E), transitions between sexual and asexual re-277 production (and vice-versa) happen at a maximum rate τ and are affected by population density. 278 We assume that individuals cannot distinguish between those that are reproducing asexually 279 and those that are potential mating partners (sexual reproduction), so population assessment is 280 based on total population size. After assessment of population density, a proportion of strategy 281 D individuals transition from D-(asexual) to D+ (sexual) and vice-versa. Similarly, a pro-282 portion of strategy E individuals transition from E-(asexual) to E+ (sexual) and vice-versa.

283
The net number of transitions from asexual (D-/E-) to asexual ( is calculated by the following equation: Equivalently, the net number of transitions from sexual (D+/E+) to asexual (D-/E-) repro-286 duction (T s− p i ,t ) is calculated by the following equation: where N s− p i ,t is the patch-specific number of individuals in strategy D/E reproducing asexually (D-288 /E-) and N s+ p i ,t is the patch-specific number of individuals in strategy D/E reproducing sexually 289 (D+/E+).

290
It is important to note that we assume no cost associated with the ability to transition in 291 strategy D/E and that these strategies are only reproducing sexually in the ancestral population 292 because population size is above the Allee saturation. We assumed τ = 0.2 across all simulations.
where µ p j→i ,t is the effective migration rate from patch p j to patch p i and µ p i→j ,t is the effective 298 migration rate from patch p i to patch p j . 299 We made migration stochastic by defining the effective migration rate according to the 300 following equation: where µ * p i→j is the mean migration rate, ρ is the magnitude of the stochasticity and N (0, 1) is 302 a standard normal deviate. We assumed µ * p S→N = 0.01, µ * p N →S = 0.001 and ρ = 0.1 across all 303 simulations.

304
Due to the stochastic environmental changes affecting population dynamics throughout time, 306 the population structure at the last time step cannot accurately describe the success of different 307 strategies in invading and establishing in the North patch. We therefore used the long-term 308 temporal mean population size of each reproductive strategy as a measurement of their ecological 309 success. Furthermore, we compared the long-term temporal mean with the initial short-term (N s p S ,0 = 0.5 · K p S ). Because the initial subpopulation size of each strategy in the ancestral 332 population is above the Allee threshold (β · K p S ), all strategies have the same initial growth 333 rate in the South patch. All simulations were performed on R version 3.5.2 [47]. considerably from that in the South patch (e.g., Figure S4). Despite fluctuations in population is much smaller than in the South patch (K p N = 0.2K p S ) and migration from the North patch to 368 the South patch occurs at a lower rate than in the opposite direction, changes in the population 369 structure in the South patch occurs at a much lower rate (e.g., Figure S5). This lower rate of is moderate (Figures S10-S13). In the long term, even when the Allee effect is strong, sexuals is low enough to overcome the Allee effect ( Figure S16). In the long term, sexuals outperform 417 parthenogens even when their sensitivity to environmental stress (φ C = 0.9) approaches that of 418 parthenogens, although in such cases the difference in ecological success is small ( Figure S17). If Note that φ represents the relative difference in sensitivity to environmental stress between apomictic (strategy A; φ A = 1.0) and automictic (strategy B; φ B ) parthenogenesis, or between parthenogenesis (strategy A; φ A = 1.0) and sexual reproduction (strategy C; φ C ).
In a deterministic simulation, with a constant environmental stress level in the marginal 3). Note that φ represents the relative difference in sensitivity to environmental stress between apomictic (strategy A; φ A = 1.0) and automictic (strategy B; φ B ) parthenogenesis, or between parthenogenesis (strategy A; φ A = 1.0) and sexual reproduction (strategy C; φ C ).
to parthenogens (φ C approaches φ A ) increases the temporal dominance of parthenogens in the 428 marginal habitat (N A /N C > 1; Figure 5). However, this effect is weaker in the long term.

429
In general, sexuals become less successful than parthenogens when β is large (e.g., β = 0.3), 430 demanding a greater population size in order to reach Allee saturation (β · K p N ), and α is Note that φ represents the relative difference in sensitivity to environmental stress between apomictic (strategy A; φ A = 1.0) and automictic (strategy B; φ B ) parthenogenesis, or between parthenogenesis (strategy A; φ A = 1.0) and sexual reproduction (strategy C; φ C ). Note that φ represents the relative difference in sensitivity to environmental stress between apomictic (strategy A; φ A = 1.0) and automictic (strategy B; φ B ) parthenogenesis, or between parthenogenesis (strategy A; φ A = 1.0) and sexual reproduction (strategy C; φ C ).
Obligate parthenogenesis vs. facultative parthenogenesis from the core habitat [20]. Geographic parthenogenesis can also result from selection for re-483 sistance and high fecundity in a facultatively parthenogenetic metapopulation occupying an 484 environmental gradient, limiting mating at the edge of the habitat and forcing females to re-485 produce asexually, generating a female bias [40]. However, the ecological success of asexual 486 reproduction certainly depends on the environmental conditions and the characteristics of the 487 competing reproductive strategy. A model exploring the frozen niche variation hypothesis [48] 488 indicated that, despite its two-fold advantage relative to sexual reproduction [2], asexual clonal 489 reproduction may have its invasion probability reduced due to a fast accumulation of delete-490 rious mutations in the initially small clonal population [49], decreasing its initial advantage of 491 local adaptation. Clonal invasion probability can be reduced even further when sexual popu-492 lations have many niche phenotypes, requiring the occurence of a beneficial mutation for niche 493 exploration within the asexual clonal population before individuals can invade a niche occupied 494 by sexuals [49]. Additionally, it has been proposed that metapopulation dynamics in marginal 495 habitats may be detrimental to sexuals because they may suffer more intensely from genetic 496 bottlenecks and inbreeding [16].

497
Although most studies focus on the coexistence or elimination of reproductive strategies as 498 a result of active selection leading to evolutionary change, it has been shown theoretically that 499 stochastic demographic processes may lead to geographic structure in the distribution of sexual 500 and asexual morphs in recently invaded areas without having to invoke adaptive differences. is low [59].

535
This positive density-dependent reproductive rate has also been explored from a theoretical proportional to the founding level of genetic variation [68]. In the clonal plant Ranunculus 565 reptans, when introduced to previously unoccupied habitats and exposed to severe stressful 566 conditions (flood and drought), populations founded by different genetic sources increased in 567 abundance relative to populations founded by genetic monocultures [69]. Interestingly, studies 568 analyzing genetic diversity as a function of distance to the core habitat (central population) 569 show that central populations have significantly higher diversity than populations located at 570 marginal habitats [70], indicating that marginal populations may be more sensitive to environ-571 mental stress and that genetic diversity may be restricted in such habitats.

572
It has been suggested that in populations where genetic diversity is very limited (e.g., clonal 573 populations) adaptation to stressful environmental conditions can be difficult when relying on 574 new mutations, leading populations to decline or go extinct [71]. This negative relationship be-tween genetic diversity and risk of extinction affects parthenogenetic populations more strongly 576 than sexual populations but this risk can be reduced when the parthenogenetic population 577 is multiclonal instead of monoclonal [48]. cial. In the Baltic Sea, asexual recruitment seems to be very common in many macrophytic 588 populations, but sexual recruitment is not completely absent [74], supporting the idea that the 589 ability to transition between sexual and asexual reproduction is beneficial. Because of all these 590 processes related to the effect of environmental stress on population growth, it is important that 591 future studies provide empirical measurements of the correlation between genetic/phenotypic 592 diversity and population sensitivity to different types of environmental stress.

594
We used a quantitative approach to explore the ecological processes that can lead to geographic  Figure S1: Logistic growth for different population growth rates. The maximum growth rate was set to r max = 1.1 in all simulations.     Figure S6: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under low environmental stress (p p N = 0.1 and v p N = 0.1), with sexual strategies exposed to a variety of Allee effect magnitudes (α and β). Axes represent the short-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5. Figure S7: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under low environmental stress (p p N = 0.1 and v p N = 0.1), with sexual strategies exposed to a variety of Allee effect magnitudes (α and β). Axes represent the long-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5. Figure S8: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under high environmental stress (p p N = 0.1 and v p N = 0.9), with sexual strategies exposed to a variety of Allee effect magnitudes (α and β). Axes represent the short-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5. Figure S9: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under high environmental stress (p p N = 0.1 and v p N = 0.9), with sexual strategies exposed to a variety of Allee effect magnitudes (α and β). Axes represent the long-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5. Figure S10: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under different probabilities of occurrence of stressful conditions (p p N ) and maximum levels of stress (v p N ), with sexual strategies exposed to a relatively weak Allee effect (α = 0.3 and β = 0.1). Axes represent the short-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5. Figure S11: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under different probabilities of occurrence of stressful conditions (p p N ) and maximum levels of stress (v p N ), with sexual strategies exposed to a relatively weak Allee effect (α = 0.3 and β = 0.1). Axes represent the long-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5. Figure S12: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under different probabilities of occurrence of stressful conditions (p p N ) and maximum levels of stress (v p N ), with sexual strategies exposed to a relatively strong Allee effect (α = 0.1 and β = 0.3). Axes represent the short-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5. Figure S13: Pairwise competition between different reproductive strategies during invasion of a marginal habitat (North patch) under different probabilities of occurrence of stressful conditions (p p N ) and maximum levels of stress (v p N ), with sexual strategies exposed to a relatively strong Allee effect (α = 0.1 and β = 0.3). Axes represent the long-term temporal mean population sizes of competing strategies. Diagonal lines indicate identical temporal mean population sizes of competing reproductive strategies and therefore equivalent ecological successes. Sensitivity values: φ A = 1.0, φ B = 0.75, φ C = 0.5.