Range expansions of sexual versus asexual organisms: Effects of reproductive assurance and migration load

It is generally considered that sexual organisms show faster evolutionary adaptation than asexual organisms because sexuals can accumulate adaptive mutations through recombination. Yet, empirical evidence often shows that the geographic range size of sexual species is narrower than that of closely related asexual species, which may seem as if asexuals can adapt to more varied environments. Two potential explanations for this apparent contradiction considered by the existing theory are reproduction assurance and migration load. Here, we consider both reproductive assurance and migration load within a single model to comparatively examine their effects on range expansions of sexuals and asexuals across an environmental gradient. The model shows that higher dispersal propensity decreases sexuals' disadvantage in reproductive assurance while increasing their disadvantage in migration load. Moreover, lower mutation rate constrains adaptation more strongly in asexuals than in sexuals. Thus, high dispersal propensity and high mutation rates promote that asexuals have wider range sizes than sexuals. Intriguingly, our model reveals that sexuals can have wider geographic range sizes than asexuals under low dispersal propensity and low mutation rates, a pattern consistent with a few exceptional empirical cases. Combining reproductive assurance and migration load provides a useful perspective to better understand the relationships between species' mating systems and their geographic ranges.

tions for this apparent contradiction considered by the existing theory are reproduction assurance and migration load. Here, we consider both reproductive assurance and migration load within a single model to comparatively examine their effects on range expansions of sexuals and asexuals across an environmental gradient. The model shows that higher dispersal propensity decreases sexuals' disadvantage in reproductive assurance while increasing their disadvantage in migration load. Moreover, lower mutation rate constrains adaptation more strongly in asexuals than in sexuals. Thus, high dispersal propensity and high mutation rates promote that asexuals have wider range sizes than sexuals. Intriguingly, our model reveals that sexuals can have wider geographic range sizes than asexuals under low dispersal propensity and low mutation rates, a pattern consistent with a few exceptional empirical cases. Combining reproductive assurance and migration load provides a useful perspective to better understand the relationships between species' mating systems and their geographic ranges.

K E Y W O R D S
gene flow, local adaptation, mating system, migration load, reproductive assurance, species' range Given the higher adaptive potential of sexually reproducing organisms than asexually reproducing organisms, it is perhaps surprising that asexual organisms tend to occupy wider geographic ranges and occur in higher altitudes and/or higher latitudes (geographic parthenogenesis, Hörandl, 2006;Tilquin & Kokko, 2016;Vandel, 1928). To explain this apparently counter-intuitive pattern, multiple non-mutually exclusive mechanisms have been proposed (Hörandl, 2006;Tilquin & Kokko, 2016). Asexual organisms are better at invading novel areas unoccupied by conspecifics because their reproduction does not require mating partners (reproductive assurance, Guzmán et al., 2012;Hörandl, 2006;Rafajlović et al., 2017).
Asexual reproduction protects marginal populations from receiving maladapted gene flow via immigrants from core habitats (migration load, Eckert, 2002;Peck et al., 1998;Piquot et al., 1996). Asexual reproduction maintains hybridization-derived populations coping better with extreme environments in marginal habitats (Bayer, 1991;Yahara, 1990). High adaptive potentials of sexual reproduction are advantageous in core habitats where antagonistic biological interactions are frequent but less important in marginal habitats (Glesener & Tilman, 1978;Scheu & Drossel, 2007;. Male mating harassment to asexual females can become weak in low-density marginal areas, promoting asexual reproduction in these habitats (Yamamichi & Koizumi, 2020). Reviewing these hypotheses, Tilquin and Kokko (2016) suggested that the enigmatic pattern of wider range sizes of asexuals than sexuals might provide key insights to understand why sexual reproduction is common despite having multiple costs relative to asexual reproduction (Lehtonen et al., 2012;Otto, 2009). Most of these hypotheses, however, have been proposed separately, with little effort having been made to examine their relative importance in explaining the broader distributions of asexuals than sexuals.
To fill this gap, we theoretically compare the conditions that determine the effects of reproductive assurance and migration load on the range expansions of sexuals and asexuals. Whereas reproductive assurance and migration load can affect the establishment of leading-edge populations in marginal habitats, the conditions required for their effects may be different. For example, dispersal can play opposite roles between these mechanisms. On one hand, dispersal between populations should tend to alleviate the disadvantage of sexual organisms in terms of reproductive assurance, because frequent immigrants from core habitats will provide more mating partners for sexual organisms in peripheral habitats. On the other hand, between-population dispersal can generate migration load, which increases the disadvantage of sexual organisms. Furthermore, the rates of mutations could affect evolutionary potentials and migration load differently between sexual and asexual reproduction.
Higher mutation rates may favour asexual over sexual reproduction because mutations might release asexual organisms from the paucity of evolutionary potentials. At the same time, high mutation rates may increase deleterious mutations (Fouqueau & Roze, 2021;Haag & Ebert, 2004;Salathé et al., 2006), which could favour either sexual or asexual reproduction depending on the fitness effects of deleterious mutations (Vanhoenacker et al., 2018).
Here, using an individual-based model, we explore the effects of dispersal and mutation rates on reproductive assurance and migration load in determining the relative geographic range sizes of sexual versus asexual organisms. The model considers that range expansions require local adaptation to a geographic gradient of environmental conditions. Although competition between sexual and asexual lineages has been an integral part of the issues on the maintenance of sex (Otto & Lenormand, 2002), we set aside the effects of their competition to focus on the effects of reproductive assurance and migration load and simulate their range expansions separately. The model will find that dispersal propensity and mutation rates affect the conditions under which reproductive assurance or migration load leads to wider geographic range sizes of asexual than sexual organisms. Intriguingly, under low dispersal propensity and low mutation rates, the model will reveal that, unlike the well-known empirical pattern, sexuals can have a wider geographic range size than asexuals.

| ME THODS
We use a spatially explicit individual-based model to simulate the range expansion dynamics of sexually and asexually reproducing organisms. To single out the effects of reproductive assurance and migration load, we assume two specific conditions. First, competition between sexual and asexual organisms has little effect on their range expansions. We, therefore, simulate their range expansion dynamics separately. Second, sexual organisms are obligately outcrossing hermaphrodites. This condition minimizes the demographic costs of sexual organisms relative to asexual organisms. Our model is applicable to plant species pairs of a sexual hermaphrodite and an asexual sister taxon. The model may also be applicable to analogous species pairs of animals, such as freshwater snails (Johnson et al., 1995) and earthworms (Mezhzherin et al., 2021).

| Model details
The range dynamics occurs along a one-dimensional stepping stone space consisting of 100 patches (indexed from j = 1 to 100). Each single patch has the common local carrying capacity (K) of 100 individuals. Each patch is characterized by environmental conditions represented by an environmental optimum. The environmental optimum E j of patch j changes along the one-dimensional line of patches.
where E 0 is set to −50 so that E j of the central patch (j = 50) is 0.
Individual organisms have a quantitative genetic trait with the trait value of individual i denoted as z i . The match between the trait value z i of an organism i and the environmental optimum E j of the patch j in which it resides defines its local performance p i as where 2 controls the degree of match necessary to achieve high performance. As explained below, we formulate fertility and viability selection such that individual fitness increases with individual performance. Thus, increasing 2 weakens fertility and viability selection. This formulation follows a standard approach to model local adaptation and range expansions (Kirkpatrick & Barton, 1997). where M is the parameter determining mate finding difficulty which expresses the minimal density to mate successfully (i.e., critical density, sensu Gerritsen, 1980), N j is the population size of patch j, and m i is a random number following the normal distribution N 1, 0.25 2 . Thus, the probability of successful mating is stochastically determined and tends to be smaller with larger M and smaller N j . When an individual acting as female is given a mating opportunity, another individual is chosen randomly from the remaining individuals of the same population to serve as male.
In the FS model, the probability that individual i ' is chosen to be a male mating partner of individual i is Assuming no genetic dominance, the trait value of an offspring between parental individuals i and i ' is defined as recombination plus Mendelian segregation and mutational random deviation: The recombination term is denoted as the mean of parents' trait values. The segregation term s is a random variable following the normal distribution with mean 0 and standard deviation 0.5. The mutation term is b • y, in which b is a random variable following a Bernoulli trial (b = 0, 1) with the mutation rate being the probability of b = 1, and y is a random variable following the normal distribution with mean 0 and standard deviation . This formulation modifies a previous approach to modelling the genetics of sexually selected traits (Doorn et al., 1998;Noest, 1997). By separating the effects of segregation and mutation, our formulation can limit the effect of Mendelian segregation only to offspring of sexual organisms while allowing the identical effect of mutation on the offspring of sexual and asexual organisms (see below). Note that an individual can mate multiple times as a male, and up to once as a female.
For asexual organisms, mating is not required for reproduction. Therefore, no mate limitation reduces the opportunity for reproduction, and all asexual individuals can produce offspring.
The number of offspring reflects the performance of parental organisms and is defined as Rp i for individual i (rounded to the nearest decimal). Offspring inherits the trait value of its parental organism with potential mutational deviation; the offspring trait value is defined as where b and y are random numbers as defined for sexual organisms.
Unlike sexual organisms, offspring are not affected by segregation.
Note that our formulation does not assume the demographic cost of sexual reproduction. A sexual individual as a female and an asexual individual produces comparable numbers of offspring when their local performances are the same. While a sexual individual as a female contributes genetically to half of its offspring, the individual can also sire the next generation as a male.
Irrespective of sexual or asexual reproduction, organisms die after reproduction (no overlapping generation), and their offspring enter the dispersal stage. The same dispersal processes apply to sexual and asexual organisms. Dispersal distance follows the Gaussian dispersal kernel. Thus, the probability that a juvenile produced in patch j disperses to patch (j + D) is where d is the standard deviation of dispersal kernel and controls dispersal propensity. When d was small, the probability of staying in the natal patch (i.e., D = 0) was high. As d increases, the probability of

| Model analysis
To explore how dispersal and mutation rates interact with reproductive assurance and migration load to determine the range sizes of sexual and asexual organisms, we analyse the model by focusing We note that this classification does not consider the possibility that the slower expanding species eventually occupies the full landscape (from j = 0 to 100) after t = 10,000, ultimately achieving the same range size as that of the faster expanding species. However, this possibility is minimal because in most cases of patterns II (sexuals wider) and III (asexuals wider), the spatial dynamics of the slower expanding species reaches equilibrium with a limited range size.

| RE SULTS
When the range sizes attained at the end of simulations are limited These patterns show that sexual populations can achieve local adaptation more efficiently than asexual populations.
Comparing the FS and VS models, local adaptation occurs more efficiently in the VS model for both sexual and asexual reproduction ( Figure 1). This is because viability selection alleviates migra- On the other hand, the mutation rate ( ) has a large impact on the range sizes of asexual organisms in the FS and VS models, with wide range sizes realized at large mutation rates (Figure 3). This is because the mutation is the sole source of evolutionary potential for local adaptation. Stabilizing selection limits range sizes, especially when the mutation rate is small but can provide sufficient evolutionary potential for a range expansion ( = 10 −4 and 10 −5 ). Unlike sexual organisms, larger dispersal propensity (d) rather facilitates range expansions as evident in the case of = 10 −5 . This is firstly because dispersal propensity does not increase migration load in asexual reproduction. Secondly, when rare mutants arise with trait values favoured in edge patches, they are more likely to be supplied to these edge patches if dispersal propensity is larger. Comparing the FS and VS models, the VS model tends to achieve larger range sizes ( Figure 3). This is because viability selection increases the opportunity of locally adapted individuals to successfully reproduce when local populations contain more juveniles (adapted and maladapted) than sustained by local carrying capacities. Not surprisingly, mate finding difficulty (M) has no influence. Extinction of asexual organisms is not observed for the studied parameter ranges.
General patterns emerge when we compare relative range sizes between sexual and asexual organisms (Figure 4). High mutation rates ( ) are the conditions most important for asexual organisms to have wider range sizes than sexual organisms (patterns III and V, see section "Model analysis" for pattern categories). Given high mutation rates ( = 10 −3 and 10 −4 ), asexuals have wider range sizes if some degree of mate limitation (M) restricts sexual reproduction.
Even without mate limitation, sufficiently strong stabilizing selection ( 2 ) leads to larger range sizes of asexuals in the FS model, which imposes a larger migration load on sexuals than in the VS model.

| DISCUSS ION
An empirical pattern has been recognized that the geographic range size of sexual organisms is generally smaller than that of closely related asexual organisms (reviewed in Hörandl, 2009;Tilquin & Kokko, 2016). Some previous models have explained this pattern by reproductive assurance giving asexual organisms an advantage in colonization (Rafajlović et al., 2017), whereas others attributed the pattern to migration load hampering the local adaptation of sexual organisms (Hu et al., 2019;Peck et al., 1998). In this paper, we have developed a simulation model that investigates how reproductive assurance and migration load together affect the range expansions of sexual and asexual organisms along an environmental gradient.
The model has shown that reproductive assurance and migration load can together result in a wider range sizes of asexual than sexual organisms, but also that their effects depend on mutation rates and dispersal propensity. In addition, the model has revealed that, unlike F I G U R E 1 Examples of range dynamics (1st and 3rd rows) and trait distribution (2nd and 4th rows) across the environmental gradient. In range dynamics, orange areas represent populations that are locally adapted (the population mean of Rp i is larger than or equal to 1), while blue areas represent populations locally maladapted (the mean of Rp i is smaller than 1). White areas represent empty patches.  (Table 1).
the well-known empirical pattern, sexual species can have wider geographic range sizes than those of asexuals.
High mutation rates are the condition important for asexual organisms to have wider range sizes than those of sexuals because otherwise asexual organisms cannot adapt to variable environmental conditions across a wide geographic range (pattern III in Figure 4). Dispersal propensity has more complex effects on relative range sizes of sexuals and asexuals. While high dispersal propensity F I G U R E 2 Effects of dispersal propensity, mutation rate, stabilizing selection and mate finding difficulty on the average range size of sexual organisms at the end of simulation. Default values are used for unvaried parameters (Table 1).

F I G U R E 3
Effects of dispersal propensity, mutation rate, stabilizing selection and mate finding difficulty on the average range size of asexual organisms at the end of simulation. Extinction is unobserved for the studied parameter ranges. Default values are used for unvaried parameters (Table 1) These results are largely consistent with previous models showing that reproductive assurance or migration load could explain wider range sizes of asexuals than sexuals (Eckert, 2002;Guzmán et al., 2012;Hörandl, 2006;Peck et al., 1998;Piquot et al., 1996;Rafajlović et al., 2017). Unlike the previous models, our model has considered both reproductive assurance and migration load within a single setting. This advantage allows us to compare the relative importance of reproductive assurance and migration load in generating the pattern of wider range sizes of asexuals. Specifically, although reproductive assurance and migration load can work together, the effect of reproductive assurance is more robust than that of migration load, because reproductive assurance can act irrespective of fertility or viability selection but migration load can be weakened under viability selection ( Figure 4). Moreover, considering that this result is derived from our model assuming hermaphrodites as sexuals, relative effects of reproductive assurance may be even stronger if sexuals are dioecious (as in most animals) because dioecious species would suffer from greater mate limitation (Gerritsen, 1980). Intriguingly, in contrast to the well-known empirical pattern, we have found that sexual species can have wider range sizes than those of asexual species (pattern II in Figure 4). In addition, sexual and asexual species can also have comparable range sizes (patterns I and IV in Figure 4). Sexual species can have wider range sizes than asexual species when mutation rates are low. Low mutation rates decrease the evolutionary potential of asexual organisms for local adaptation, while sexually reproducing organisms can nonetheless maintain their evolutionary potential. In our model, sexual organisms can maintain evolutionary potential despite low mutation rates because sexual reproduction can integrate rare mutant trait values into the genetic variation of populations; without sexual reproduction, rare mutant trait values do not contribute to population's evolutionary adaptation unless mutant lineages become the majority of the population (Burt, 2000;Fisher, 1930;Muller, 1932;Neher et al., 2010). Given low mutation rates, wider range sizes of sexual than asexual species occur if mate limitation and migration load are weak enough. These conditions are more likely with viability selection and under high dispersal propensity because viability selection alleviating migration load mitigates the negative effect of high dispersal propensity in increasing migration load. Although our model considers only mutations affecting the trait determining local adaptation, high mutation rates may generally increase mutation load more on asexual organisms (Vanhoenacker et al., 2018). In such a F I G U R E 4 Effects of dispersal propensity, mutation rate, stabilizing selection and mate finding difficulty on the relative range sizes between sexual and asexual organisms at the end of simulation. Asexual extinction and both extinction (patterns VI and VII) are unobserved for the studied parameter ranges. Default values are used for unvaried parameters (Table 1).
case, sexuals might have larger range sizes than asexuals even when mutation rates are high.
Empirical evidence supports our prediction that high mutation rates are necessary for asexuals to have wider range sizes than sexuals ( Figure 4). In experimental and natural bacteria, asexually reproducing populations evolved to have high mutation rates (Sniegowski et al., 2000). Moreover, ploidy has been considered as a main mechanism underlying wider range sizes of asexual insects than those of closely related sexuals, because multiploidy increases the total mutation rate across the whole genome (Lundmark & Saura, 2006).
These examples indicate that asexual organisms maintaining high mutation rates satisfy the condition necessary for asexuals to have wider range sizes than sexuals.
There are also some empirical examples supporting our prediction that sexuals can have wider range sizes than asexuals. Examining seven pairs of outcrossing and selfing sister species from Collinsia, Grant and Kalisz (2020) found the outcrosser of one pair having a wider niche breadth than the selfer (Grant & Kalisz, 2020). In another example, a sexual plant species (Cortaderia selloana) is more invasive and has a wider distribution than a closely related asexual agamospermous species (Cortaderia jubata) in its introduced range (Lambrinos, 2001). Since our model prediction is that this pattern occurs if mate limitation and migration load are weak, it would be interesting to explore the causes of the corresponding empirical patterns by examining strengths of mate limitation and migration load in these systems.
A recent review on the causes of range limitation found little evidence that gene flow restricts species ranges (Kottler et al., 2021).
The review suggested that rather than migration load, the lack of gene flow into isolated populations at range edges makes them suffer from inbreeding depression and thereby limits species ranges.
This view agrees with our results that lower dispersal propensity, although weakening migration load, tightens mate limitation to restrict the range expansion of sexual organisms (Figure 2) (Case & Taper, 2000;Encinas-Viso et al., 2020). The absence of competition might be a valid assumption for comparing the range expansion of sexual and asexual species on a large spatial scale, on which abiotic conditions rather than local competition can be important for delimiting species distributions (Connell, 1980;Connor & Simberloff, 1979;Strong Jr et al., 1979). Still, if sexuals and asexuals do compete locally, outcomes of local competition can affect their distributions on larger spatial scales. We might expect that local competition can strengthen or weaken the predictions of our model. For example, pre-emptive competition would permit a species that firstly occupies a local habitat patch, facilitating the species expanding more quickly to achieve a wider range size (Kilsdonk & De Meester, 2021;Vanoverbeke et al., 2016). Alternatively, the speed of range expansion would be unimportant if over a long term, a late arriving species competitively excludes an early arriving species from local patches. Possibly, late arriving sexuals might catch up early arriving asexuals (as in Cuellar, 1977;Rafajlović et al., 2017;Soreng & Van Devender, 1989) and exclude the early occupants after antagonistic interspecific interactions (e.g. competition, disease and predation) become more important than mate limitation or migration load as envisaged in the Red Queen hypothesis (Hamilton, 1980;Hamilton et al., 1990;Van Valen, 1973) and the tangled-bank hypothesis (Bell, 1982;Ghiselin, 1974; In real systems, environments can change more abruptly between neighbouring habitat patches (Liebherr, 1988;Sandoval, 1994). Since large environmental differences between neighbouring patches can make migration load more intense, abrupt environmental changes can further decelerate the range expansions of sexuals (Bridle et al., 2019).
In conclusion, our model has shown that mutation rates, dispersal propensity and fertility versus viability selection can affect reproductive assurance and migration load to determine relative advantages of sexuals and asexuals in range expansions along an environmental gradient. In addition to clarifying the conditions for asexuals to have wider geographic range sizes than sexuals, our model has revealed the conditions under which sexuals can be better at range expansions than asexuals. These conditions, such as low mutation rates and high dispersal propensity coupled with viability selection, might also help understand why sexual reproduction is common despite having multiple costs relative to asexual reproduction (Lehtonen et al., 2012;Otto, 2009). Our results may also be useful for predicting how different mating systems affect the responses of species distributions to future environmental changes.

ACK N OWLED G M ENTS
We thank Deciding Editor Dr. Luke Holman and two anonymous reviewers for encouraging and constructive comments that greatly improved the manuscript. Open access publishing facilitated by The University of Tokyo, as part of the Wiley -The University of Tokyo agreement.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have no conflict of interest to declare.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.14161.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in "Dryad" at https://doi.org/10.5061/dryad.59zw3 r2c8, reference number dryad.59zw3r2c8.