Charline Smadi

Hélène Leman

Violaine Llaurens
email: violaine.llaurens@mnhn.fr

Looking for the right mate in diploid species: how does genetic dominance affect the spatial differentiation of a sexual trait?

Keywords: Mate preference, Heterozygote, Spatial segregation, Migration, Uneven population size

come

Looking for the right mate in diploid species: how dominance relationships affect population differentiation in sexual trait?

Introduction

Understanding processes leading to biological diversification is a central question in evolutionary biology. Traits may diverge neutrally because of geographic barriers limiting gene flow, or simple isolation by distance due to limited dispersal, resulting in genetic and phenotypic differentiation [START_REF] Lande | Genetic variation and phenotypic evolution during allopatric speciation[END_REF]; [START_REF] Slatkin | Gene flow and the geographic structure of natural populations[END_REF]). Mating preferences may also drive divergence in targeted traits, with assortative mating promoting local fixation of the most abundant phenotype. Migration behaviours can then be affected by such mating preferences, because mate searching can stimulate dispersal when the preferred trait is locally rare. Altogether, mate preference has been shown to have a great impact on both local polymorphism [START_REF] Payne | Sexual selection, space, and speciation[END_REF] and spatial structure (M'Gonigle and FitzJohn (2010)) for targeted traits. Nevertheless the genetic architecture of the trait under sexual selection may also influence the evolution of such population differentiation. Linkage disequilibrium between loci controlling the adaptive trait and the preference trait is known to favour divergence via a Fisherian run-away process [START_REF] Ronald | The genetical theory of natural selection: a complete variorum edition[END_REF]). Dominance relationships among differentiated alleles may also greatly influence the spatial distribution of different phenotypes [START_REF] John R Pannell | haldane's sieve' in a metapopulation: sifting through plant reproductive polymorphisms[END_REF]). The effective migration of advantageous alleles is favoured when they are dominant: immigrant alleles entering a new population will mostly occur in heterozygous state so that they will be picked up by positive selection only if they are expressed. Recessive adaptive alleles, scarcely expressed, are more likely to be lost by genetic drift because of their neutrality when in a heterozygous state. This Haldane's sieve effect [START_REF] Burdon | A mathematical theory of natural and artificial selection, part v: selection and mutation[END_REF]) predicts a greater effective migration of dominant adaptive alleles as compared to recessive ones. For example, in the polymorphic locus of sporophytic self-incompatibility where rare alleles benefit from increased reproductive success, migration of dominant alleles has been shown to be more effective than migration of recessive ones [START_REF] Mikkel | Evolutionary dynamics of sporophytic self-incompatibility alleles in plants[END_REF]). Dominance among alleles may thus play an important role in the dynamics of spatial dif-ferentiation of traits under sexual selection. Here we investigate the influence of dominance on the spatial differentiation of a trait determining both mating and migration behaviours. Our model is an extension of a previously described model assuming haploid individuals [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]). Two populations, linked by migration, are assumed. A mating trait, encoded by one locus has two consequences: (1) encountering pairs mate more often when both individuals display the same phenotype for this trait and (2) the migration rate of an individual is proportional to the frequency of the other phenotype in its patch: individuals are more prone to move if they have difficulties finding a suitable mate in their patch. This hypothesis is relevant for all organisms with active mate searching (e.g. patrolling behaviours in butterflies [START_REF] Jugovic | Movement, demography and behaviour of a highly mobile species: A case study of the black-veined white, aporia crataegi (lepidoptera: Pieridae)[END_REF]) or mate sampling in bowerbirds [START_REF] Albert | Complex mate searching in the satin bowerbird ptilonorhynchus violaceus[END_REF])), but also in organisms where gametes are involved in dispersal (e.g. sea urchins (Crimaldi (2012)) or plants [START_REF] Millar | Extensive long-distance pollen dispersal and highly outcrossed mating in historically small and disjunct populations of acacia woodmaniorum (fabaceae), a rare banded iron formation endemic[END_REF])), and may travel large distances before encountering a suitable mate. This original hypothesis linking mating preference and migration contrasts with assumptions found in classical models of speciation (see [START_REF] Gavrilets | Models of speciation: Where are we now[END_REF]) for a review) where the preference traits are generally not directly linked to migration behaviour. In the haploid model studied in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]), this preference behaviour can lead to spatial differentiation of the trait between the two populations with fixation of different phenotypes in the two populations despite migration, and without any local adaptation. More precisely, if, initially, one phenotype is predominant in one patch and the other phenotype predominant in the other patch, then spatial differentiation occurs regardless of the migration rate. Given the importance of dominance on the migration of alleles, here we extend this model to a diploid case with dominance relationships between alleles, and explore how dominance may modulate the dynamics of spatial differentiation on this mating trait.

Materials and methods

In this section we provide a general description of the model of population dynamics. We also describe the mathematical methods used to derive analytical results as well as the conditions used for numerical simulations. More details can be found in the Supplementary Materials.

We consider a population of hermaphroditic diploid individuals characterized by (1) a single phenotype controlled by their genotype at one bi-allelic locus (A and a), and (2) their occupancy of one of two patches (1 and 2). The number of individuals AA, Aa and aa in the two patches follow a dynamical system, which can be obtained as a limit of a stochastic multitype birth and death process with competition in continuous time ( [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]) and see Appendix A). In particular, population sizes in the patches are varying and generations are overlapping.

The phenotype of individuals influences both (1) their mating and (2) migration behaviours. (1) Individuals having the same phenotype have a higher probability of mating (see [START_REF] Jiang | Assortative mating in animals[END_REF]) for a recent review of the mechanisms of assortative matings in animals). ( 2) This mate preference also influ-ences migration from one patch to another: individuals carrying a phenotype at low frequency within patch have a greater migration because we assumed migration to be promoted by the local lack of suitable mates. Therefore, the migration rate depends on the frequency of individuals carrying a different phenotype within their patch. Examples of animals migrating to find suitable mates are well documented [START_REF] Schwagmeyer | Scramble-competition polygyny in an asocial mammal: Male mobility and mating success[END_REF]; [START_REF] Höner | Female mate-choice drives the evolution of male-biased dispersal in a social mammal[END_REF]. A migration mechanism similar to the one presented here has been studied in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]) and in [START_REF] Payne | Sexual selection, space, and speciation[END_REF]) in a continuous space model.

Five parameters are needed to describe the dynamics:

• b is the minimal individual birth rate. It corresponds to the rate at which an individual gives birth if there is no individual with the same phenotype in its patch.

• β ≥ 1 is the sexual preference. Individuals encounter each other randomly, and two individuals with similar phenotypes have a higher probability to mate (see Section 2.1 for details).

• p describes the individuals' ability to migrate. The effective migration rate of an individual is the product of p and of the proportion of 'not preferred' individuals (see Section 2.2 and Appendix A for details).

• d is the individual death rate in the absence of competition.

• c represents the competition for food or space exerted by an individual on another one of the same patch. An additional individual death rate results from the total competition exerted on this individual.

There is no geographically variable selection: the parameters b, p, d, and c do not depend on the patch, thus the two patches have the same ecological characteristics. However, depending on the patch's composition, the individual's behaviour in terms of reproduction and migration may vary in time and space.

Dominance

We assume Mendelian segregation so that any parent transmits each of its alleles with probability 1/2. To study the effect of dominance on population differentiation, we contrast two opposite scenarios: complete co-dominance and complete dominance. Dominance observed in natural populations can differ from these extreme cases, however by studying the limits of the dominance continuum we cover the possible population dynamics. In the complete dominance scenario, individuals with genotypes AA and Aa have the same phenotype, A, whereas individuals with genotype aa have the phenotype a. In the co-dominance scenario, heterozygotes express an intermediate phenotype between the two homozygotes and we thus consider two possibilities. (1) Preference expressed towards heterozygotes will be intermediate between assortative and disassortative mating (COD 1) (β for pairs (AA, AA), (aa, aa) and (Aa, Aa), (β + 1)/2 for pairs (AA, Aa) and (aa, Aa), and 1 for pairs (AA, aa)) and migration rate varies accordingly because the decision to leave the patch depends on the lack of preferred partners in the patch. (2) Heterozygotes express no preference towards any partners (COD2), the preference parameter is thus β for pairs (AA, AA) and (aa, aa), and 1 for all other pairs. Because of this lack of preference, heterozygotes have no reason to discriminate against the individuals in their patch and thus do not migrate.

Model

Our model can be precisely described as follows. We denote the number of individuals of genotype g in the patch i at time t for any g ∈ {AA, Aa, aa}, i ∈ {1, 2} and t ≥ 0 by z g , i(t). Moreover we denote the total number of individuals in the patch i at time t by

N i (t) := z AA,i (t) + z Aa,i (t) + z aa,i (t).
The parameter β measures the strength of the sexual preference: β = 1 means no preference, and a large β indicates a strong preference. Since we are only investigating assortative mating here, β will always be larger than one. We denote by p β (g, g ) the preference between two individuals with genotypes g and g , respectively. They differ according to the model considered (see Table 1) and always belong to [1, β]. Finally, the parameter p describes the maximum migration rate of an individual. The migration rate of an individual is the product of this parameter p, and the proportion of unpreferred mates in its patch, and the function of preferences p β (., .) between genotypes.

The equations governing the evolution of the population were obtained as follows. Individuals are assumed hermaphroditic and diploid. At a rate B > bβ, they will reproduce as female and look for a mate. They will encounter uniformly an individual reproducing as male from the population. The probability that they actually mate and produce an offspring is bp β (g, g )/B, where the functions p β are detailed in Table 1.

In the patch 1, AA individuals reproducing as female will generate:

• an offspring of type AA in the patch 1 at a rate

Bz AA,1 z AA,1 + z Aa,1 + z aa,1 bβ B z AA,1 + bp β (AA, Aa) B z Aa,1 2 .
Indeed, as we consider a Mendelian segregation, parents of genotypes AA and Aa generate an offspring of type AA with a probability 1/2 and an offspring of type Aa with a probability 1/2.

• an offspring of type Aa in the patch 1 at a rate

Bz AA,1 z AA,1 + z Aa,1 + z aa,1 b B z aa,1 + bp β (AA, Aa) B z Aa,1 2 .
Since individuals are hermaphrodites and can reproduce through both male and female pathways, we finally get the following equation for the population dynamics on a patch i ∈ {1, 2} (where j is the complement of i in {1, 2}):

                                                                                                 żAA,i = b N i βz 2 AA,i + p β (AA, Aa)z AA,i z Aa,i + p β (Aa, Aa) 4 z 2 Aa,i -(d + cN i )z AA,i + p g∈G p β (AA, AA) -p β (AA, g) β -1 z g, j N j z AA, j - z g,i N i z AA,i żAa,i = b N i p β (Aa, Aa) 2 z 2 Aa,i + p β (AA, Aa)z Aa,i z AA,i + p β (aa, Aa)z Aa,i z aa,i + 2z AA,i z aa,i -(d + cN i )z Aa,i + p g∈G p β (Aa, Aa) -p β (Aa, g) β -1 z g, j N j z Aa, j - z g,i N i z Aa,i żaa,i = b N i βz 2 aa,i + p β (aa, Aa)z aa,i z Aa,i + p β (Aa, Aa) 4 z 2 Aa,i -(d + cN i )z aa,i + p g∈G p β (aa, aa) -p β (aa, g) β -1 z g, j N j z aa, j - z g,i N i z aa,i
(1) where G := {AA, Aa, aa}. Note that individual migration rates are always between 0 and p. The equations under the three dominance hypotheses are then obtained by replacing the preference functions p β by their value, as summarized in Table 1 (see Appendix A for a presentation of the full equations in the different models). (COD1) (COD2) (DOM) p β (AA, AA)

β β β p β (AA, Aa) (β + 1)/2 1 β p β (Aa, Aa) β 1 β p β (aa, Aa) (β + 1)/2 1 1 p β (aa, aa) β β β p β (AA, aa) 1 1 1

Mathematical analysis

We carry out a mathematical analysis of the dynamical systems governing the population evolution, fully described in Supplementary Materials. Here, we only show some fixed points of the systems, their stability, and obtain some convergence results in the case without migration. We used the theory of dynamical systems (in particular Lyapunov functions and the Local Center Manifold Theorem) as well as the theory of polynomial functions.

Simulations

To illustrate the evolution of the population sizes under our model, we used the software Mathematica Inc.. We also performed numerical simulations of the dynamical systems presented in Appendix A. We investigated the solution of the dynamical systems for b = 2, d = 1, c = 0.5 and different values of the migration rate p and of the preference coefficient β. For each value of β and p, we solved the dynamical systems for 10 000 different initial conditions chosen as follows. First, we set the size of the population in patch 1: we considered 100 values of sizes regularly distributed between 1 and 2 * (2 * bd)/c. Then, we set the size of the population in patch 2 such that the difference between the two sizes was 0.01. This allowed us to reduce the number of parameters explored without changing the results. Finally, for each pair of sizes, we examined 100 initial conditions randomly chosen. For each patch i, we set uniformly at random the proportion p A,i of allele A using a uniform random variable between 0.5 and 1 in the case where the majority allele was A or between 0 and 0.5 in the other case. The proportion of Aa-individuals in patch i was fixed randomly using a uniform random variable between 0 and 2 * min(p A,i , 1p A,i ). The proportion of AA-individuals and aa-individuals in each patch can be easily deduced.

For each initial condition, we numerically solved the three dynamical systems using a finite difference method. We used a discretization time step h equals to 0.005. We assumed that a stationary state was reached as soon as the norm of the difference of the solution between two time steps was lower than ε = 10 -6 . Different values of h and ε were tested and the chosen pair of values provided the best trade-off between algorithm rapidity and result accuracy. Once the stationary state was found, we considered that the final population in a patch was monomorphic when the proportion of one of the alleles was larger than 99%. Otherwise, it was considered to be a polymorphic equilibrium.

Results

To highlight the influence of dominance on the spatial dynamics of the trait, we first contrast our findings with results obtained in a previous study [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]) where individuals were haploids.

Haploid model

In the haploid model, preference functions were p β (A, A) = p β (a, a) = β and p β (A, a) = 1. The system may converge to two different types of equilibria depending on initial conditions. Equilibria can be expressed using ζ, the equilibrium population size in a patch when there is only one type of individual (A or a),

ζ := βb -d c . (2)
Let us denote by z α,i (t), α ∈ {A, a}, i ∈ {1, 2}, t ≥ 0 the α population size in the patch i at time t. Then

• If z A,i (0) > z a,i (0), i = 1, 2, the population sizes converge to (z A,1 , z a,1 , z A,2 , z a,2 ) = (ζ, 0, ζ, 0). • If z A,1 (0) > z a,1 (0) and z A,2 (0) < z a,2 (0) the population sizes converge to (z A,1 , z a,1 , z A,2 , z a,2 ) = (ζ, 0, 0, ζ).
The same conclusions hold by symmetry if we replace A by a in the first line and 1 by 2 in the second one. Hence when conditions are initially asymmetrical (more A individuals in one patch and more a individuals in the other patch) and when p is positive, regardless of its value, it is enough to entail the end of gene flow between the two populations and to generate two differentiated populations.

Diploid models with co-dominance

We then compared the results obtained in haploid populations with the diploid model assuming co-dominance between alleles, i.e. when heterozygotes express an intermediate phenotype between either homozygotes. This assumption exhibits high similarity with the haploid case, although the behaviour of heterozygous individuals displaying intermediate phenotypes might influence model outputs. We contrasted two hypotheses:

(1) mating success between homozygotes and heterozygotes was half that of matings between identical genotypes (COD1), (2) heterozygotes had no preferences and were not preferred by homozygotes (COD2) (preference parameter 1 for any encountering pair with an Aa individual). We then investigated the equilibrium reached in both populations using a mathematical analysis and simulations assuming different preference coefficients (β) and migration rates (p).

System without migration

When there is no migration (p = 0) we are able to give necessary and sufficient conditions on the initial numbers of individuals with genotypes AA, Aa and aa in both patches for the system to converge to the different fixed points, in both codominant models. Recall that z αα ,i (t) denotes the number of individuals with genotype αα (AA, Aa or aa) in the patch i (1 or 2) at time t. As there is no migration, it is enough to consider the patch 1. Under hypothesis (COD1) the system follows the equations:

                             żAA = b N βz 2 AA + β + 1 2 z AA z Aa + β 4 z 2 Aa -(d + cN)z AA żAa = b N β 2 z 2 Aa + β + 1 2 z Aa (z AA + z aa ) + 2z AA z aa -(d + cN)z Aa żaa = b N βz 2 aa + β + 1 2 z aa z Aa + β 4 z 2 Aa -(d + cN)z aa ,
where N = z AA + z Aa + z aa . This system admits two stable fixed points, (z

AA,1 = ζ, z Aa,1 = 0, z aa,1 = 0) (fixation of al- lele A) and (z AA,1 = 0, z Aa,1 = 0, z aa,1 = ζ) (fixation of allele a)
, where we recall that ζ has been defined in (2), and one unstable fixed point with persistence of all three genotypes (see Appendix C.1). If z AA,1 (0) > z aa,1 (0), allele A becomes fixed, and the numbers of individuals converge to the stable equilibrium (ζ, 0, 0). If z AA,1 (0) < z aa,1 (0), the exact same conclusion holds with a replacing A. Finally, if z AA,1 (0) = z aa,1 (0), the system converges to the unstable equilibrium. Under hypothesis (COD2), the equations driving the population dynamics are:

                           żAA = b N βz 2 AA + z AA z Aa + 1 4 z 2 Aa -(d + cN)z AA żAa = b N z 2 Aa + z Aa (z AA + z aa ) + 2z AA z aa -(d + cN)z Aa żaa = b N βz 2 aa + z aa z Aa + 1 4 z 2 Aa -(d + cN)z aa .
We get the same result, except that the unstable fixed point is different (see Appendix C.2 for details).

Fixed points of the system with migration

When there is a migration between the two patches (p > 0), the dynamics are much more complex due to the increase in dimensionality and we were unable to obtain convergence results analytically. However we were able to describe some of the fixed points and determine their stability. There are four fixed points with monomorphic populations in both patches: fixation of A in both patches (z AA,1 = ζ, z Aa,1 = 0, z aa,1 = 0, z AA,2 = ζ, z Aa,2 = 0, z aa,2 = 0), fixation of a in both patches (0, 0, ζ, 0, 0, ζ), or fixation of different alleles in the two patches, (ζ, 0, 0, 0, 0, ζ) or (0, 0, ζ, ζ, 0, 0). For the two codominant models, the first two fixed points are stable for all parameters values, and the two last fixed points are stable if p < p crit = bβ(β -1)/2 and unstable if p > p crit . The fact that p crit is proportional to b is expected, as we could reduce the number of parameters in Equation ( 1) by taking b = b/p, d = d/p, c = c/p and p = 1. p crit is increasing with β. This is also expected as a higher β means a higher local advantage due to sexual preference for the major allele in a patch. Hence a higher p is needed to counteract this advantage by enhancing gene flow between patches. The form of p crit as a function of β however is difficult to infer, except the fact that it becomes null when β = 1. Indeed, in this last case, the model becomes completely neutral and an infinity of equilibria are possible (total population size (bd)/c in each patch and any proportions of AA, Aa, and aa individuals). This result contrasts with the haploid case, where the fixed points with a genotype in each patch were stable for all values of the migration parameter p (see [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]). This may be explained by the fact that the migration of heterozygotes has a major impact on the population behaviour, as demonstrated below using numerical simulations. In the figure displaying the behaviour of the model in the codominant case, we have indicated the curve p crit = bβ(β -1)/2 (see Fig. 1).

Conditions for differentiated populations

Under both assumptions regarding co-dominance, differentiated populations with fixation of different phenotypes in the two populations can emerge only when the frequencies of allele a are asymmetrical at initial state, with frequency of allele a smaller than 1/2 in one population and larger than 1/2 in the other population (see Appendix E.2). As shown in Appendix D.1 and Appendix D.2, this equilibrium cannot be reached when p > bβ(β -1)/2. However, an initial asymmetry and the condition p < bβ(β -1)/2 are not enough to ensure differentiated populations, as can be seen from Fig. 1.

The number of simulations exhibiting differentiation increases when the preference coefficient β increases as expected, and decreases when the migration rate p increases. Indeed when migration increases or preference decreases, this equilibrium becomes unstable (see black curve in Fig. 1) and is no longer reached. This contrasts with the dynamics of the model when individuals are haploid, where differentiated populations emerge as soon as the frequencies of allele a are asymmetrical at initial state, regardless of the value of the migration strength, Figure 1: Conditions for the emergence of differentiated populations with fixation of allele A in population 1 and a in population 2, assuming hypothesis (1) regarding codominance (outcomes are very similar under the alternative hypothesis of codominance). Initial conditions are asymmetric, with more than 50% of A allele in patch 1, and more than 50% of a allele in patch 2. Color indicates the percentage of simulations where A fixes in population 1 and a in population 2. Black line shows the limit of the stability of the equilibrium with differentiated populations, p crit = bβ(β -1)/2. p, and preference β > 1. However, the negative effect of migration on population differentiation is observed under both assumptions regarding co-dominance, including hypothesis (2) where heterozygotes never migrate. This indicates that their presence may be enough to maintain the migration of both A and a homozygotes across populations when p is large, even if they do not move themselves.

3.2.4. Fixation of a single phenotype throughout both populations Under hypothesis (COD1), when migration increases, the two populations tend to be more homogeneous, leading to the fixation of a single allele throughout both populations, whatever the initial frequencies. The identity of the fixed allele depends on the initial frequency because of positive frequencydependent selection triggered by homogamy. In cases where initial frequencies of allele a are asymmetrical in the two populations, we indeed observed fixation of allele a in 50% of simulations when migration increases (Fig 2a-b). As the system is symmetrical in a and A, in the other half of simulations, the fixation of allele A was observed (data not shown). When allele a initially predominates in both patches, the fixation of allele a is observed in all simulations for hypothesis (COD1) regarding heterozygote behaviour (Fig. 2c). When heterozygotes express a preference for themselves and migrate, their migration leads to a fast equalization of the numbers of individuals with the same genotype in the two populations (see Fig. 3) (i.e. n AA,1 (resp. n Aa,1 , n aa,1 ) very close to n AA,2 (resp. n Aa,1 , n aa,2 )). Once this equalization is reached, the migration does not influence the dynamics since the numbers of emigrants and immigrants of each population are the same. Both patches evolve as if they were isolated. From the study of the system without migration in Section 3.2.1, we know that there is no polymorphic stable equilibrium in one isolated patch. As a consequence, the same allele gets fixed in both patches.

Assuming no preference and no migration of heterozygotes (hypothesis (COD2)) leads to the overall fixation of one of the two alleles only when migration is limited or preference is high (Fig. 2d). Otherwise in many simulations, polymorphism persists in both populations.

Polymorphic equilibria

The dynamics leading to these equilibria are highly nonmonotonic and migration persists at equilibrium, which explains the difficulty of studying the model analytically. The presence of individuals of type Aa maintains migration of AA and aa individuals at equilibrium (see Fig. 4 and5). According to numerical simulations, stable polymorphic equilibria are of the following form:

-one 'large' population (with a size larger than (βbd)/c) with homozygote type in large majority, AA or aa respectively, -one 'small' population (with a size smaller than (bd)/c) with a large proportion of individuals of type aa (resp. AA) and Aa.

This result is rather unexpected since (βbd)/c is the maximal equilibrium population size of an isolated population whereas (bd)/c is smaller than the minimal one. Indeed, the maximal birth rate of any individual is βb, which is its birth rate when surrounded by individuals with the same phenotype, whereas its minimal birth rate is b if surrounded by individuals with different phenotypes. In all cases, its death rate is the sum of its baseline death rate, d, and of the competition death rate which

(0) = z AA,2 (0) = z Aa,2 (0) = 0.1, z Aa,1 (0) = 2, z aa,1 (0) = 1 and z aa,2 (0) = 2.
is equal to the product of c and of the total population size in its patch. As a consequence, in an isolated patch without any migration, this leads to a maximal (resp. Moreover, equilibria with polymorphism are observed only when the migration rate p is high with respect to the preference coefficient β (see Fig. 2d), which reinforces the idea that polymorphism is maintained by a trade-off between migration and selection.

Dominance between alleles

We now assume total dominance of the allele A, so that Aa heterozygotes display the same phenotype and behaviour as AA homozygotes.

System without migration

When there is no migration (p = 0), the system admits two stable fixed points, (z AA,1 = ζ, z Aa,1 = 0, z aa,1 = 0) (fixation of allele A) and (z AA,1 = 0, z Aa,1 = 0, z aa,1 = ζ) (fixation of allele a), and one unstable fixed point with persistence of all three genotypes (see Appendix C.3). In this case we were only able to give a sufficient condition on the initial number of individuals of different types for the system to converge to the stable fixed point characterizing the fixation of allele A: if z AA,1 (0) ≥ z aa,1 (0), the solution converges to the stable equilibrium (ζ, 0, 0).

Fixed points of the system with migration

When there is a migration between the two patches (p > 0), once again the dynamics is much more complex and we were unable to obtain convergence results. However we were able to describe some of the fixed points and determine their stability (see Appendix D.3). There are four fixed points with monomorphic populations in both patches, which are the same as in the case of codominance: fixation of A in both patches (z AA,1 = ζ, z Aa,1 = 0, z aa,1 = 0, z AA,2 = ζ, z Aa,2 = 0, z aa,2 = 0), fixation of a in both patches (0, 0, ζ, 0, 0, ζ), or fixation of different alleles in the two patches, (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0). The first two fixed points are stable for all the parameters values. The two last fixed points admit five negative eigenvalues and one null eigenvalue and we were not able to assess their stability. Since the recessive allele a has no influence on the phenotypes of Aa heterozygotes, it behaves neutrally with respect to mating and migration behaviour when occurring in heterozygotes. This neutral behaviour may substantially complexify the two-populations dynamics observed here.

Numerical simulations were then performed to study the influence of preference and migration parameters on the equilibria reached.

Fixation of the dominant allele A

In many simulations assuming an asymmetrical initial state (more a alleles in one patch and more A alleles in the other patch), the fixation of the dominant allele A was observed, when Figure 5: Proportions of different genotypes and migration rates between the two patches at equilibrium. z AA,1 ∼ 0.06, z Aa,1 ∼ 0.64, z aa,1 ∼ 0.96, z AA,2 ∼ 0.45, z Aa,2 ∼ 0.36, z aa,2 ∼ 2.19. Total population in patch 1, 1.77, total population in patch 2, 3. Same parameters, initial conditions and colors as in Figure 4 preference was strictly larger than 1 (Fig. 6a). Dominance of allele A makes the phenotype A (displayed by both AA and Aa genotypes) more frequent therefore promoting its fixation through assortative mating advantage. When allele A was initially the minor allele however, its fixation happened in many fewer simulations, because the large initial frequency of genotypes aa displaying phenotype a overrides the dominance effect (Fig. 6b). Interestingly, although increasing migration promotes the fixation of the dominant allele A through its homogenizing effect, increasing preference tends to weaken this fixation. This non-trivial effect may stem from the limitation of migration when a homogeneous population emerges: individuals matching the predominant phenotype within a patch migrate less. Therefore, when preference is strong and populations are initially differentiated, equilibria with fixation of different alleles in the two populations might be frequent, despite the frequencydependent advantage of the phenotype carried by the dominant allele.

Conditions for differentiated populations

As observed in Fig. 7a, simulations with initial differentiation mainly result in the fixation of different alleles in populations 1 and 2 when the preference coefficient β increases and migration strength is not too high. However, it is worth noting that differentiated populations persist for values of the migration parameter p much higher than in both codominant cases (see Fig. 1). When allele a is initially the major allele in both patches, simulations mostly lead to the fixation of the recessive allele a (data not shown). Note that around 25% of simulations lead to differentiated populations whereas in the haploid version of the model studied in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]) differentiated populations could emerge only when the initial state was asymmetrical. Note also that in some parameter regions (small preference parameter β or high migration strength p) no simulation leads to differentiated populations. This may indicate that equilibria (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0) (see Section 3.3.2) are unstable for such parameters.

Polymorphic equilibria

In some cases, even if initial conditions are asymmetrical (more a alleles in one patch and more A alleles in the other patch), genetic polymorphism can persist in one or both populations. However, the populations are almost phenotypically monomorphic in both patches: almost only aa individuals in one patch, and almost only AA and Aa individuals in the other one (see Fig. 8 and9). As a consequence, individuals reproduce at their maximal birth rate βb and do not migrate, and the population sizes in both patches are close to their carrying capacity (bβd)/c.

Constant migration

In order to understand better the role of frequency dependent migration, we explored the behaviours of the different genotypes when the migration rate is constant, equal to p and identical for all individuals.

In this case, there are three possible types of equilibria z (see Appendix F for the proof): (1) Either z = 0

(2) Or all the coordinates of z are positive (the three genotypes are present in the two patches)

(3) Or there is only one type of individual, AA or aa, and the population size is the same in the two patches,

ζ := βb -d c .
Although a strict division of the population into two monomorphic subpopulations with different allele (A in a patch and a in the other one) is not possible, there are cases where AA is in large majority (more than 90%) in one patch, and aa in the other patch. This behavior can be observed under both codominant hypotheses (simulations not shown) as well as under the dominant hypothesis, as presented in Fig. 10. However, note that the preference parameter (β) has to be very large and the migration parameter (p) very small. The ratio between the preference parameter and the migration parameter has to be much larger than in Section 3.3.4 with our initial model.

Actually, this condition on the ratio between the preference parameter and the migration rate is necessary to observe equilibria which are not of type (3) described above. In Fig. 11, we show results assuming strict dominance between alleles. Hence, the particular behaviours of these models with reduced migration are observed only for large values of β (β larger than 3 for p = 0.5 for example), which contrasts substantially with our initial models. The particular migration behaviour assumed in our initial model, which depends on the local composition of potential mates, indeed substantially promote population differentiation.

Discussion

Effect of ploidy on population differentiation

Altogether, these results show that simulations of diploid organisms assuming codominance between alleles depart from equilibria observed in the haploid model. Notably, migration limits the probability of population differentiation in the diploid model which was not observed in haploids. Population differentiation is highly sensitive to initial conditions, because alleles with high frequencies are strongly advantaged by assortative preferences, favouring fixation of the predominant phenotype within populations. Because migration depends on the number of individuals with different phenotypes, population differentiation is quickly achieved in the haploid model because of a rapid decrease of migration as soon as differentiation starts, until a complete isolation of the two populations, which fixed different alleles. In the diploid model however, when the two alleles are codominant, the presence of heterozygotes with intermediate phenotypes promotes migration even when populations are initially differentiated, and frequently leads to the fixation of a single allele throughout both connected populations. The discrepancy observed here between haploid and diploid assumptions highlights the need to consider the effect of ploidy in spatiallystructured models of trait evolution, because the presence of an intermediate phenotype can interfere in the differentiation process.

Effect of intermediate phenotypes' behaviour on polymorphism

Depending on the assumption regarding codominance, the resulting equilibrium slightly differed. The absence of preference of heterozygotes for their own phenotypes (hypothesis (COD2)) can promote polymorphism of the trait under sexual selection, notably when populations are initially uneven. These heterozygotes do not migrate and reproduce equally with any genotypes. When they initially occur in significant proportion within populations, they promote migration of homozygous genotypes and limit fixation of the initially predominant allele. This contrasts with the co-dominance drawn from hypothesis (COD1), where heterozygotes can migrate and are half as preferred by homozygotes. In this case, fixation of a single allele throughout both populations is always achieved when the same allele is predominant in both populations initially.

Heterozygote behaviour is therefore a key parameter in the dynamics of population differentiation in mating traits. However, empirical data on mating behaviour of intermediate phenotypes are scarce. In Heliconius butterflies, wing colour pattern is known to be an important visual cue for mate choice. In the specific case of Heliconius heurippa, which displays a red and white colour pattern, which can be obtained by crossing its red sister species H. melpomene to the white sister species H. cydno, assortative mate preferences have been demonstrated [START_REF] Mavárez | Speciation by hybridization in heliconius butterflies[END_REF]). This assortative mating has been hypothesized to favour the emergence of the species H. heurippa, putatively created by hybridization events between H. melpomene and H. cydno [START_REF] Jiggins | Hybrid trait speciation and heliconius butterflies[END_REF]). Similarly, the golden-crowned manakin species (Lepidothrix vilasboasi) has recently been demonstrated to have emerged from an hybrid speciation between the snow-capped (Lepidothrix nattereri) and opal-crowned (Lepidothrix iris) manakins of the Amazon basin, leading to an intermediate phenotype (Barrera-Guzmán et al. ( 2017)). Female choice has been suggested to play a key role in the persistence of this new phenotype, driving the emergence of a new species.

We therefore hope that our predictions on the important impact of co-dominant heterozygote behaviour on the evolution of mating trait differentiation will motivate further empirical research on mating and migration behaviour of intermediate phenotypes.

Evolutionary consequences of dominance for population differentiation and speciation

The fixation of the recessive allele a when initially predominant in both populations is rarer when the alternative allele is dominant as compared to codominant. Dominant alleles spread among populations more easily than codominant ones because of a Haldane's sieve effect [START_REF] Burdon | A mathematical theory of natural and artificial selection, part v: selection and mutation[END_REF]). Assuming strict dominance, fixation of the dominant allele throughout both populations is thus frequent, and favoured when predominant in both populations (data not shown).

However, the invasion of the dominant haplotype throughout both populations is limited when initial populations display uneven proportions of the two alleles, leading to either (1) population differentiation or (2) persistent polymorphism within populations.

(1) Dominance may reinforce population differentiation, because heterozygotes display the dominant phenotype and therefore rapidly increase the number of individuals carrying this phenotype within one population, therefore causing a decrease in migration between populations. Consequently, population differentiation is more frequently observed when one allele is dominant as compared to co-dominant, even when assuming a high migration rate.

(2) Depending on the initial distribution of heterozygotes among populations, polymorphism can also be maintained within a population where the dominant phenotype is frequent: heterozygotes then display the preferred phenotype and therefore do not suffer from mate rejection and scarcely migrate. This last result is in accordance with a recent paper [START_REF] David M Schneider | Diploid versus haploid models of neutral speciation[END_REF], where the authors explore the consequences of dominance at loci involved in genetic incompatibilities on the dynamics of speciation in a spatially-explicit individually centered model. They observed that the distance between mates necessary for a spatial mosaic of species to emerge needed to be more restricted in model assuming diploidy with strict dominance as compared to haploid model. This highlights how dominance may modulate spatial differentiation and the emergence of well-separated species. Altogether this highlights the need for diploid models of speciation, and should stimulate empirical comparisons of speciation dynamics driven by adaptive or sexually-selected traits displaying contrasting dominance relationships.

Aknowledgements

This work was funded by the Chair "Modélisation Mathématique et Biodiversité" of VEOLIA-Ecole Polytechnique-MNHN-F.X to HL and CS, the Young Researcher ANR DOMEVOL (ANR-13-JSV7-0003-01) and the Emergence program of Paris city council to VL.

Appendix A. Definition of the models

In these Supplementary Materials, we will provide more details of the models, address the question of fitnesses, and prove the results presented in the main text.

Recall that we denote the number of individuals of type α in the patch i at time t for any α ∈ {AA, Aa, aa}, i ∈ {1, 2} and t ≥ 0 by z α,i (t). Moreover we denote the total number of individuals in the patch i at time t by

N i (t) := z AA,i (t) + z Aa,i (t) + z aa,i (t).
Finally, we recall that the parameter

ζ := bβ -d c
is the equilibrium size of a monomorphic AA or aa population.

It will be a characteristic quantity in many equilibria for the three dynamical systems.

Recall that the models have been described in Section 2.2. In the codominant case, we consider two different models, depending on the heterozygotes behaviour:

(1) Either the preference expressed by an individual is an 'average preference' of its alleles (β for pairs (AA, AA), (aa, aa) and (Aa, Aa), (β + 1)/2 for pairs (AA, Aa) and (aa, Aa), and 1 for pairs (AA, aa)) and the migration rate follows the same rule,

(2) Or heterozygotes express no preference and thus do not migrate (in this case only the homozygotes express a preference towards individuals of the same genotype). The preference parameter is thus β for pairs (AA, AA) and (aa, aa), and 1 for the other pairs.

For the first model (COD1, we obtain: The second codominant model yields:

                                                                 żAA,i = b N i βz 2 AA,i + β + 1 2 z AA,i z Aa,i + β 4 z 2 Aa,i -(d + cN i )z AA,i -p z aa,i + z Aa,i /2 N i z AA,i + p z aa, j + z Aa, j /2 N j z AA, j żAa,i = b N i β 2 z 2 Aa,i + β + 1 2 z Aa,i (z AA,i + z aa,i ) + 2z AA,i z aa,i -(d + cN i )z
                                                         żAA,i = b N i        βz 2 AA,i + z 2 Aa,i 4 + z Aa,i z AA,i        -(d + cN i )z AA,i -p z aa,i + z Aa,i N i z AA,i + p z aa, j + z Aa, j N j z AA, j żAa,i = b N i        (z AA,i + z aa,i )z Aa,i + z 2 Aa,i 2 + 2z AA,i z aa,i        -(d + cN i )z Aa,i żaa,i = b N i        βz 2 aa,i + z 2 Aa,i 4 + z Aa,i z aa,i        -(d + cN i )z aa,i
p z AA,i + z Aa,i N i z aa,i + p z AA, j + z Aa, j N j z aa, j .

(A.2) In the dominant case, only aa individuals express the phenotype a. Other individuals are of phenotype A. The population dynamics in this case are:

                                                             żAA,i = bβ N i z AA,i + 1 2 z Aa,i 2 -(d + cN i )z AA,i -p z aa,i N i z AA,i + p z aa, j N j z AA, j żAa,i = b N i (βz Aa,i + 2z aa,i ) z AA,i + z Aa,i 2 -(d + cN i )z Aa,i -p z aa,i N i z Aa,i + p z aa, j N j z Aa, j żaa,i = b N i βz 2 aa,i + z aa,i z Aa,i + β 4 z 2 Aa,i -(d + cN i )z aa,i
p z AA,i + z Aa,i N i z aa,i + p z AA, j + z Aa, j N j z aa, j , (A.3) Notice that the dynamical systems governing the population dynamics in the co-dominant and the dominant cases can be obtained as large population limits of stochastic individual based models (see [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]).

To be precise, let K be a large parameter that gives the order size of the population. The microscopic population is represented by the process (Z K α,i (t)) α∈{AA,Aa,aa},i∈{1,2} where Z K α,i (t) gives the size, divided by K, of the α-population in patch i at time t. The reproduction and migration mechanisms are described in Section 2.2 and the death rate of any individual in patch i is given by

d + c Z K AA,i (t) + Z K Aa,i (t) + Z K aa,i (t) .
Under the assumption that the sequence of initial conditions (Z K α,i (0)) α∈{AA,Aa,aa},i∈{1,2} K∈N converges (in probability) when K goes to infinity, the sequence of stochastic functions (Z K α,i (t)) α∈{AA,Aa,aa},i∈{1,2} , t ∈ [0, T ] K∈N also converges (in probability for the uniform convergence) to the trajectory of the deterministic models presented above. In Fig. A.12 and A.13, we present a realisation of the trajectory of Z K aa,1 for different values of K and under the parameters values presented in Fig. 3 and4. Notice that K = 10000 is enough for the limiting deterministic model to be a very good approximation of the stochastic process. The mechanism of mating preference presented in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]) is similar to the ones classically used in the ecology literature (see for instance [START_REF] Gavrilets | On the evolution of premating isolation after a founder event[END_REF]; [START_REF] Matessi | Long-term buildup of reproductive isolation promoted by disruptive selection: how far does it go?[END_REF]; [START_REF] Bürger | Intraspecific competitive divergence and convergence under assortative mating[END_REF]; Servedio (2010)) and references therein). The unusual form of the equations comes from the fact that we model a varying size population evolving in continuous time and with overlapping generations, whereas classical models consider discrete non overlapping generations models with an infinite population size (see Section 'Discussion of the model' in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]) for a detailed comparison). As in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF] and in order to compare both continuous and discrete models, we give the probabilities that the individuals with genotype g mate with any individual with genotype g in the deme i at time t: p β (g, g )z 2 g,i N i p β (g, AA)z AA,i + p β (g, Aa)z Aa,i + p β (g, aa)z aa,i .

Finally, to ensure the survival of the population and that the total population size remains bounded, we make the following assumptions:

b > d > 0 and c > 0.

We recall that we are interested in the case of assortative mating, which means that β ≥ 1.

Proof. To begin with, we describe the different stable fixed points in the codominant case. We recall that the dynamics of the system in one patch is given by the system of equations:

                             żAA = b N βz 2 AA + β + 1 2 z AA z Aa + β 4 z 2 Aa -(d + cN)z AA żAa = b N β 2 z 2 Aa + β + 1 2 z Aa (z AA + z aa ) + 2z AA z aa -(d + cN)z Aa żaa = b N βz 2 aa + β + 1 2 z aa z Aa + β 4 z 2 Aa -(d + cN)z aa , (C.1) where N = z AA + z Aa + z aa .
The fixed points (z AA , z Aa , z aa ) are solutions to the following system of equations:

                             b N βz 2 AA + β + 1 2 z AA z Aa + β 4 z 2 Aa = (d + cN)z AA b N β 2 z 2 Aa + β + 1 2 z Aa (z AA + z aa ) + 2z AA,1 z aa = (d + cN)z Aa b N βz 2 aa + β + 1 2 z aa z Aa + β 4 z 2 Aa = (d + cN)z aa . (C.2)

Monomorphic equilibria

We first check that if z αα = 0 for α ∈ {A, a}, then z Aa is necessarily equal to 0. Hence we obtain the two following monomorphic fixed points: They are all negative under our assumptions. If z Aa > 0 then from the first and the last equations in (C.2) we see that necessarily z AA > 0 and z aa > 0. Hence we look for a fixed point with the three coordinates positive. From the first and the last equations in (C.2) we find that

(d + cN)z AA z aa = z aa (βz 2 AA + β + 1 2 z AA z Aa + β 4 z 2 Aa ) = z AA (βz 2 aa + β + 1 2 z AA z Aa + β 4 z 2
+ β + 1 2 z aa z Aa + β 4 z 2 Aa = 1 z Aa       β 2 z 2 Aa + β + 1 2 z Aa       z 2 Aa 4z aa + z aa       + z 2 Aa 2       . This implies β -1 2       z aa + z 2 Aa 4z aa       = 0,
and contradicts the fact that all the fixed point coordinates are positive.

Case z aa = z AA

The fixed point is thus of the form (x, y, x), with x > 0 and y > 0. Equalizing the first equation of (C.2) divided by x and the second one divided by y, we find

βx + β + 1 2 y + β 4 y 2 x = β 2 y + (β + 1)x + 2 x 2 y ⇔ x 3 y 3 + 1 2 x 2 y 2 - 1 4 x y - β 8 = 0.
Let us set δ := x y > 0. The polynomial function P

(X) = X 3 + 1 2 X 2 -1 4 X -β
8 , has only one positive root and this latter belongs to ]1/2, +∞[. Indeed by taking the first derivative of P we can check that P(X) is increasing until x = -1/2, then decreasing until X = 1/6 and then increasing again. As P(-1/2) = (1 -β)/8 < 0, we conclude that there is only one positive root. Finally we can check that P(1/2) < 0. Then using the fact that x = δy and the second equation of (C.2), we conclude that there exists only one equilibrium with positive coordinates which is (δξ, ξ, δξ) where ξ :

= b(δ + β/2) -d (2δ + 1)c ,
and δ is the unique positive solution of P(X) = 0. Using again that P(δ) = 0, we can write the eigenvalues of the Jacobian Matrix at (δξ, ξ, δξ) as

bβ 4δ (2δ -1), λ 2 , λ 3 ,
where λ 2 and λ 3 are complex numbers. As δ ∈]1/2, +∞[, this ensures that the fixed point (δξ, ξ, δξ) is unstable. We now prove that its stable manifold is of dimension 2 and is exactly the set

A = = {z ∈ (R * + ) 3 , z AA = z aa }.
By subtracting the third equation of (C.1) to the first one, we get

żAA -żaa = (z AA -z aa ) bβ -d -cN -b β -1 2 z Aa N , (C.4)
and we deduce that the set A = is invariant under the flow defined by (C.1). Thus, let us assume that the initial condition belongs to A = . Hence for any t ≥ 0, z AA (t) = z aa (t). Moreover,

d dt z AA (t) δz Aa(t) -1 = - 2b δz 2 Aa N       z 3 AA + z Aa 2 z 2 AA - z 2 Aa 4 z AA - βz 3 Aa 8       .
Since P(δ) = 0, δz Aa is a root of the polynomial function with respect to the variable z AA of the right hand side. Thus

d dt z AA (t) δz Aa(t) -1 = - 2b z Aa N z AA δz Aa -1       z 2 AA + 1 2 + δ z Aa z AA + βz 2 Aa 8δ       . (C.5)
The polynomial function of degree two of the r.h.s in (C.5) is non-negative. This ensures that the function

W 1 (z) = ln z AA δz Aa - 1
is a Lyapounov function for the dynamical system (C.1) restricted on A = . Theorem 1 in LaSalle (1960) implies that the flow z(t) converges to {(δξ, ξ, δξ)}, the largest invariant set included in {z AA = δz Aa } ∩ A = . This implies that A = is included in the stable manifold of (δξ, ξ, δξ).

Let us now deal with the solution outside A = . First note that by (C.4), the two sets

A > = {z ∈ (R * + ) 3 , z AA > z aa } and A < = {z ∈ (R * + ) 3 , z AA < z aa }
are two invariant sets under the dynamical system (C.1). Since the system (C.1) is symetric with respect to z AA and z aa , we only have to deal with one of the previous sets, the dynamics in the other one being symmetric. In what follows, we study the dynamics in the set A > . On the set A > , the Lyapunov function W 2 (z) = ln z AA + z aa + βz Aa z AAz aa is well defined and is clearly non-negative. Moreover, its derivative is where N = z AA + z Aa + z aa .

In this section, we will prove the following result:

Lemma 3. The dynamical system (C.6) admits exactly three fixed points:

• Two stable fixed points (ζ, 0, 0) and (0, 0, ζ).

• One unstable fixed point b(β + 1) -2d 4c

       β + 1 -1 β + 1 + 1 , 2 β + 1 + 1 , 1        .
• (0, 0, ζ, 0, 0, ζ) (equilibrium with genotype aa in both patches) is stable and there exists a neighborhood V aa of (0, 0, ζ, 0, 0, ζ) such that any solution starting from V aa converges exponentially fast to the equilibrium,

• (ζ, 0, 0, ζ, 0, 0) (equilibrium with genotype AA in both patches) is stable and there exists a neighborhood V AA of (ζ, 0, 0, ζ, 0, 0) such that any solution starting from V AA converges to the equilibrium with a rate t → 1 t , • (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0) are two equilibria whose Jacobian Matrices admit five negative eigenvalues and a null eigenvalue.

In the case of the two last equilibria, we were not able to prove theoretically their stability. Simulations of the solution seem to show that the stability depends on the values of the parameters.

Proof. Equilibrium (0, 0, ζ, 0, 0, ζ):

The fixed point (0, 0, ζ, 0, 0, ζ) is the easiest to study, as all the eigenvalues of the Jacobian Matrix at this point are negative: Hence, the equilibrium (0, 0, ζ, 0, 0, ζ) is stable.

Equilibrium (ζ, 0, 0, ζ, 0, 0):

The second fixed point, (ζ, 0, 0, ζ, 0, 0), is more involved to study, since the Jacobian Matrix at this points admits four negative eigenvalues, and 0 as an eigenvalue with multiplicity two: (0, 0, -bβ, -bβ + d, -bβ + d, -bβ -2p).

As a consequence, we have to go further to get the stability of this fixed point. The Jacobian Matrix at this fixed point is

                        -bβ + d -bβ + d -2bβ + d -p 0 0 p 0 0 2b 0 0 0 0 0 -bβ -p 0 0 p 0 0 p -bβ + d -bβ + d -2bβ + d -p 0 0 0 0 0 2b 0 0 p 0 0 -bβ -p                         , (D.1)
and its eigensystem is (0; (0, 0, 0, -1, 1, 0)) (0; (-1, 1, 0, 0, 0, 0)) Hence if we make a translation such that the fixed point has coordinates (0, 0, 0, 0, 0, 0), and introduce a new basis via the following matrix:

                              0 -1 --2bβ 2 +2bβ+dβ-2d βd 0 1 -2β 2 b 2 -2βb 2 -βdb+2db+6βpb+4p 2 -2d p (bβ+2p)(d+2p) 0 1 -2 β 0 0 2b bβ+2p 0 0 1 0 0 -1 -1 0 --2bβ 2 +2bβ+dβ-2d βd 1 0 2b(d-bβ) (bβ+2p)(d+2p) --2bβ+d-2p d+2p 1 0 -2 β 0 0 -2b
We notice that the coefficients of P are positive. Hence any real root is negative. This proves that all the roots are negative in the case where the three roots are real numbers. The second possibility is that P admits one negative real root (denoted by -ρ) and two complex conjugate roots (denoted by φe iθ and φe -iθ , φ ∈ R + , θ ∈ [0, 2π)). We aim at showing that cos θ < 0, which will ensure that all the roots of P have negative real parts.

With these notations, P can be rewritten as follows P(X) = X 3 + (ρ -2φ cos θ)X 2 + φ(φ -2ρ cos θ)X + ρφ 2 .

As the coefficients of P are positive, we get ρ > 2φ cos θ and φ > 2ρ cos θ, which yields 2 cos θ < ρ φ < 1 2 cos θ , and thus 4 cos 2 θ < 1.

But we also have the following series of inequalities: 0 < b 3 β(2+7(β-1)β)+b 2 (2+7β(3β-2))p+4b(5β-1)p 2 +6p 3 = AE -2F = (ρ -2φ cos θ)φ(φ -2ρ cos θ) -2ρφ 2 = φ 2 ρ(4 cos 2 θ -1) -2φ cos θ(φ 2 + ρ 2 ) ≤ -2φ cos θ(φ 2 + ρ 2 ).

This implies that cos θ < 0, as expected. Unfortunately, as we do not know explicitely the roots of P, we cannot apply Theorem 2.12.1 in [START_REF] Perko | Differential equations and dynamical systems[END_REF] to get the behaviour of the solutions of (C.6) in the center manifold and deduce the stability of the fixed point (0, 0, ζ, ζ, 0, 0).

Appendix E. Some properties of the system with migration

Figure 2 :

 2 Figure 2: Conditions of fixation of allele a in both populations. Columns represent simulations assuming hypotheses (1) and (2) regarding heterozygote behaviour. Rows differ in initial conditions. First row: Asymmetrical frequency of allele a (more frequent in population 2); Second row: Frequency of allele a greater than 0.5 in both populations. The colors indicate the percentage of simulations where a fixes in both populations.

Figure 3 :

 3 Figure 3: Number of individuals in both patches under hypothesis (1). β = 1.1 and p = 5. Colors: the dynamics of the number of individuals in patch 1 (resp.2) are represented in red (resp. black), the dynamics of the number of individuals with genotype AA (resp Aa, aa) are represented using a solid line (resp. dashed line, dotted line) The initial conditions are z AA,1 (0) = z AA,2 (0) = z Aa,2 (0) = 0.1, z Aa,1 (0) = 2, z aa,1 (0) = 1 and z aa,2 (0) = 2.

  minimal) equilibrium equal to (βbd)/c (resp. larger than (bd)/c). Actually, polymorphism can be maintained through an equilibrium between growth and migration. Indeed, denoting by N the population size in the less populated patch, we deduce that the number of births per time unit in this patch always stays above b * N, whereas the number of deaths is (d + c * N) * N. As N is smaller than (b-d)/c, the number of deaths can not exceed (d + (bd))N = b * N and there are more births than deaths in the less populated patch. Since the equilibrium is maintained, there is a continuous excess flux of migration towards the other patch. Only AA and aa individuals migrate, which explains why the proportion of Aa individuals remains relatively high. On the contrary, as the population size of the most populated patch is larger than (βbd)/c, there are more deaths than births in it, but the latter constantly receives individuals of types AA and aa from the less populated patch, which maintains the polymorphism.

Figure 4 :

 4 Figure 4: Number of individuals carrying each genotype and their migration rates under hypothesis (2). β = 1.1 and p = 5, with the parameters, (bβ-d)/c = 2.4. The initial conditions are z AA,1(0) = z AA,2 (0) = z Aa,2 (0) = 0.1, z Aa,1 (0) = 2,z aa,1 (0) = 1 and z aa,2 (0) = 2. Colors: (a) the dynamics in patch 1 (resp. 2) are represented in red (resp. black), the dynamics of the number of individuals with genotype AA (resp Aa, aa) are represented using solid (resp. dashed, dotted) lines; (b) the migration from patch 1 to patch 2 (resp. 2 to 1) is drawn in red (resp. black), the migration of aa-individuals (resp. AA) is represented using dotted (resp. solid) lines.

Figure 6 :

 6 Figure 6: Conditions of fixation of the dominant allele A in both populations. Columns differ in initial conditions. First column: Asymmetrical frequency of allele a (more frequent in population 2); Second column: Frequency of allele a greater than 0.5 in both populations. The color indicates the percentage of simulations where A fixes in both populations.

Figure 7 :

 7 Figure 7: Conditions for the emergence of differentiated populations (with fixation of allele A in population 1 and a in population 2). Columns differ in initial conditions. First column: Asymmetrical frequency of allele a (more frequent in population 2); Second column: Frequency of allele a greater than 0.5 in both populations. The colors indicate the percentage of simulations where A fixes in population 1 and a in population 2.

Figure 8 :

 8 Figure 8: Population sizes and migration rates under the dominance hypothesis. β = 1.5 and p = 5. The initial conditions are z aa,1 (0) = 0.5 and 0.1 for the initial number of individuals with the other genotype. Colors: (a) the dynamics in patch 1 (resp.2) are represented in red (resp. black), the dynamics of the number of individuals with genotype AA (resp Aa, aa) are represented using solid (resp. dashed, dotted) lines; (b) the migration from patch 1 to patch 2 (resp. 2 to 1) is drawn in red (resp. black), the migration of aa-individuals (resp. AA) is represented using dotted (resp. solid) lines.

Figure 9 :

 9 Figure 9: Proportions of different genotypes in the two patches at equilibrium. z AA,1 ∼ 0.01, z Aa,1 ∼ 0.009, z aa,1 ∼ 3.96, z AA,2 ∼ 3.45, z Aa,2 ∼ 0.51, z aa,2 ∼ 0.02. Total population sizes in patch 1 and 2 very close to (bβd)/c = 4. Same parameters, initial conditions and colors as in Figure 8

Figure 10 :

 10 Figure 10:Conditions for the emergence of differentiated populations. Columns differ in initial conditions. First column: Asymmetrical frequency of allele a (more frequent in population 2); Second column: Frequency of allele a greater than 0.5 in both populations. The colors indicate the percentage of simulations where the patch 1 is filled with more than 90% of allele A in population 1 and the patch 2 with more than 90% of allele a. Note that the scales for β and p are really different from those of Fig.7

Figure 11 :

 11 Figure 11: Conditions for equilibria with all coordinates positive. Columns differ in initial conditions. First column: Asymmetrical frequency of allele a (more frequent in population 2); Second column: Frequency of allele a greater than 0.5 in both populations. The colors indicate the percentage of simulations with equilibria which are not of type (3).

Figure A. 12 :

 12 Figure A.12: Comparison of the probabilistic and the deterministic COD1 models. Colored curves: realisations of the trajectory (Z K aa,1 (t), t ≥ 0) for five different values of K given in the legend; black curve: deterministic trajectory. The parameters are the ones of Fig. 3.

Figure A. 13 :

 13 Figure A.13: Comparison of the probabilistic and the deterministic COD2 models. Colored curves: realisations of the trajectory (Z K aa,1 (t), t ≥ 0) for five different values of K given in the legend; black curve: deterministic trajectory. The parameters are the ones of Fig. 4.

  1), -bβ + d .

  AA + z aa + βz Aa ) [βz Aa (z AA + z aa )] ≤ 0, and this derivative is equal to 0 on the set O > = {z ∈ A > , z Aa = 0}. The largest invariant set with respect to (C.1) including in O > is {z ∈ A > , z Aa = z aa = 0} and it is obvious that any trajectory starting from this invariant set converges to (ζ, 0, 0). Using Theorem 1 in[START_REF] Joseph P Lasalle | Some extensions of liapunov's second method. Circuit Theory[END_REF], it is sufficient to conclude that any trajectory starting from A > converges to (ζ, 0, 0).-(d + cN)z AA żAa = b N (βz Aa + 2z aa ) z AA + z Aa 2 -(d + cN)z Aa

(

  -bβ, bbβ, -bβ + d, -bβ + d, -bβ -2p, bbβ -2p).

Table 1

 1

	: Preference functions assuming dominance hypothesis (COD1),
	(COD2) and (DOM).

  Aa or z aa = z AA .

				z 2 Aa /4z aa , we get
				1 z aa	βz 2 aa
				Aa ),
	which yields
	(z AA -z aa )(z AA z aa -	z 2 Aa 4	) = 0.	(C.3)
	Hence, either 4z AA z aa = z 2
	Case 4z AA z aa = z 2 Aa
	By expressing in two different ways (d + cN)N thanks to
	the second and last equations in (C.2) and replacing z AA by

  bβ -2p; -(2b 2 β(β -1) + bd(2 -β) + 6bβp -2d p + 4p 2 ) (d + 2p)(bβ + 2p)

						,	2b bβ + 2p	,
	-1,	2bβ -d + 2p d + 2p	+	2b(d -bβ) (bβ + 2p)(d + 2p)	, -	2b bβ + 2p	, 1 .
									-bβ;	2bβ(β -1) + 2d -βd βd	, -	2 β	, 1,	2bβ(β -1) + 2d -βd βd	, -	2 β	, 1

(-bβ + d; (0, 0, 0, 1, 0, 0)) (-bβ + d; (1, 0, 0, 0, 0, 0)) -

Appendix B. On fitness proxies

The question of defining a fitness proxy for diploid individuals with a density dependent sexual reproduction is tricky. Fitnesses are mostly defined in the case of haploid individuals with clonal reproduction. It thus essentially consists in computing the exponential growth rate of a certain type of individual (for instance mutants) in a resident population. It can also be computed for individuals in a population with more than two types. Hence it may be necessary to take into account the current population sizes.

A first solution in our case could be to compute the relative rates at which AA and aa individuals take part in an event of reproduction, as an AA (resp. aa) individual necessarily transmits an allele A (resp. a). However, such a fitness proxy would not allow us to take properly into account the role of the heterozygote individuals, especially when one of the alleles is dominant.

Very few fitness proxies to the diploid case have been available so far (but see [START_REF] Roze | Inbreeding depression and the evolution of dispersal rates: a multilocus model[END_REF]; [START_REF] Ravigné | Selective interactions between short-distance pollen and seed dispersal in self-compatible species[END_REF]; [START_REF] Parvinen | A novel fitness proxy in structured locally finite metapopulations with diploid genetics, with an application to dispersal evolution[END_REF]; [START_REF] Parvinen | On fitness in metapopulations that are both size-and stage-structured[END_REF]). These papers essentially deal with a meta-population and dispersers. We will borrow some ideas of [START_REF] Parvinen | A novel fitness proxy in structured locally finite metapopulations with diploid genetics, with an application to dispersal evolution[END_REF]) and of branching process theory (see [START_REF] Athreya | Branching processes[END_REF] for instance) to propose a fitness proxy in our case.

In [START_REF] Parvinen | A novel fitness proxy in structured locally finite metapopulations with diploid genetics, with an application to dispersal evolution[END_REF]), the authors devise a fitness proxy for a set of metapopulation models defined in continuous time. It corresponds to the expected number of mutant dispersers produced by a local mutant population initiated by one mutant disperser and is computed via the principal eigenvalue of a matrix taking into account the fact that a mutant allele may be in a homozygote or an heterozygote mutant.

Our case is more involved however, because the birth rates are type and density dependent. Our strategy is thus the following: we consider an environment set by the current population sizes in a patch, construct a multitype branching process with this environment, and define relative fitnesses as the long term proportions of individuals. It has to be understood that such fitnesses are relative. We aim at finding which genotype is locally favored in a patch among AA, Aa, and aa, and not giving a quantitative fitness value to each genotype.

To construct this fitness proxy, let us first recall that an interpretation of the birth rates is that half of the time individuals reproduce as a female (they choose their mate according to their preference), and half of the time as a male (they are chosen by the female). For a fixed value of population sizes in a patch (z AA , z Aa , z aa ) we may thus compute the rate at which an individual with a given genotype σ gives birth to an individual of a given genotype τ, b στ .

In our model, these rates have the following expressions (where z αα , (α, α ) ∈ {A, a} 2 is the αα -population size in the patch considered, ᾱ is the complementary of α in {A, a}, and

These rates can be seen as the birth rates of a continuous time multitype branching process. If we introduce the matrix

We know that this matrix has a positive maximal eigenvalue λ, a left and a right eigenvectors u and v associated to λ with positive coordinates such that u.v = 1 and u.1 = 1.

Then the total population size of the branching process grows exponentially with a Malthusian parameter λ and the proportions of the populations (z AA , z Aa , z aa ) converge to u. These proportions could be a proxy for the fitnesses of the different genotypes in a patch at a given moment.

Appendix C. Behaviour of the system without migration

Let us first study the behaviour of the system when there is no migration. In this case the two patches have independent dynamics, and it is enough to study one patch.

Appendix C.1. First codominant case

We will prove the following result: Lemma 1. If p = 0, there are two stable and one unstable fixed points in the patch i for i ∈ {1, 2}: (ζ, 0, 0), (0, 0, ζ) and (δξ, ξ, δξ), where δ is the unique positive root of the polynomial functional

Moreover, we have the following asymptotic behaviours for the dynamical system (A.1):

• If z AA,i (0) > z aa,i (0), then the solution converges to the stable fixed point (ζ, 0, 0)

0), then the solution converges to the unstable fixed point (δξ, ξ, δξ).

• If z AA,i (0) < z aa,i (0), then the solution converges to the stable fixed point (0, 0, ζ)

Appendix C.2. Second codominant case Using the same ideas and the same proofs as in the previous section Appendix C.1, we prove the following lemma for the codominant model when heterozygotes are not preferred and do not migrate: Lemma 2. If p = 0, there are two stable and one unstable fixed points in the patch i for i ∈ {1, 2}:

where δ is the unique positive root of the polynomial functional

Moreover, we have the following asymptotic behaviours for the dynamical system (A.1):

, then the solution converges to the stable fixed point (ζ, 0, 0)

, then the solution converges to the unstable fixed point (δ ξ , ξ , δ ξ ).

• If z AA,i (0) < z aa,i (0), then the solution converges to the stable fixed point (0, 0, ζ)

The proof is really similar except that we use the Lyapunov function

to end it.

Appendix C.3. Dominant case

We recall that the dynamics of the system in one patch is given by the system of equations:

Proof. The fixed points (z AA , z Aa , z aa ) are solutions to the following system of equations:

Again we can show that if z αα = 0 for α ∈ {A, a}, then necessarily z Aa = 0, and we find the two following monomorphic equilibria (ζ, 0, 0) and (0, 0, ζ)

Hence the equilibrium (0, 0, ζ) is stable. To find the stability of the equilibrium (ζ, 0, 0), we will use a Lyapunov function.

First notice that the derivative of the total population size satisfies:

and the derivative of the difference z AAz aa satisfies:

In particular, the set A := {z AA ≥ z aa } is stable. Indeed, if we consider a point where z AA = z aa then the previous derivative is positive and we stay in the set A. Moreover in this set

We deduce thatln((z AA -z aa )/N) is a Lyapunov function which cancels out of the set {z aa = 0} ∪ {z AA = z Aa = 0}. Beside the only fixed point in this set is (ζ, 0, 0). Applying Theorem 1 in LaSalle (1960), we deduce that any solution to (C.6) starting from A converges to the fixed point (ζ, 0, 0).

Let us now find and study the positive equilibrium. By subtracting the second equality to the first one in (C.7), we get

which gives the equality

Hence, by subtracting the third to the first inequality in (C.7), we find

Finally, by using (C.8), we deduce

which leads to z AA + z Aa = z aa .

(C.9)

From (C.8) and (C.9), we get

If we inject these inequalities in the derivative of z Aa , we obtain

.

We deduce that the positive equilibrium we are looking for, if it exists, has necessarily the following coordinates

Conversely, we can check that this point is indeed an equilibrium. The eigenvalues of the Jacobian matrix at this point are (2β( β + 1 + 1) 2 ) -1 times the roots of the polynomial

We deduce that there are two negative and one positive eigenvalues. This equilibrium is thus unstable. As a conclusion, in the dominant case, the only stable fixed points are (ζ, 0, 0) and (0, 0, ζ).

Appendix D. Stability of the monomorphic fixed points in the system with migration

As the systems (A.1) and (A.3) are complex, we are not able to derive an analytical form for all the fixed points. As a consequence in the next two subsections we focus on the fixed points which are monomorphic in one patch, (ζ, 0, 0, ζ, 0, 0), (0, 0, ζ, 0, 0, ζ), (ζ, 0, 0, 0, 0, ζ), and (0, 0, ζ, ζ, 0, 0).

Appendix D.1. First codominant case

We will prove that for the first codominant model the monomorphic fixed points have the following properties: Lemma 4. There are four fixed points which are monomorphic in every patch: , 0, 0, ζ, 0, 0) and(0, 0, ζ, 0, 0, ζ) are stable. • (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0) are stable if 2p < bβ(β-1) and unstable if 2p > bβ(β -1).

Since the alleles a and A are codominant, the two fixed points (ζ, 0, 0, ζ, 0, 0) and (0, 0, ζ, 0, 0, ζ) have the same characteristics. The eigenvalues of the Jacobian Matrix at these fixed points are

They are all negative, and we conclude that (ζ, 0, 0, ζ, 0, 0) and (0, 0, ζ, 0, 0, ζ) are two stable fixed points.

In the codominant case, the two fixed points (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0) have the same properties since the alleles a and A are codominant. The eigenvalues of the Jacobian Matrix at these fixed points are

All the eigenvalues, except the fourth one, are negative. The fourth eigenvalue is negative only if the migration parameter p is small enough:

Hence the eigenvalues (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0) are stable only if 2p < bβ(β -1).

Appendix D.2. Second codominant case

Lemma 4 is still valid for the second codominant case, although the eigenvalues of the fixed points are not the same:

for the fixed points (ζ, 0, 0, ζ, 0, 0) and (0, 0, ζ, 0, 0, ζ), and

for the fixed points (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0).

Appendix D.3. Dominant case

In the dominant case, the two fixed points (ζ, 0, 0, ζ, 0, 0) and (0, 0, ζ, 0, 0, ζ) are not symmetrical, as the two alleles A and a do not have the same role. We prove the following Lemma.

Lemma 5. There are four fixed points which are monomorphic in every patch:    .

From the first equation in (D.3) we know that in a neighbour-hood of (0, 0), Ẋ = G 1 (X, h(X)). Hence

But we notice that the solutions ( x1 , x2 ) of

Moreover, using that

Hence x 1 + x 2 has a constant sign and converges to 0 with a rate t → 1 t in a neighbourhood of the fixed point (0, 0, 0, 0, 0, 0). But by inversion of the matrix (D.2), we see that

This ends the proof of the stability of the fixed point (ζ, 0, 0, ζ, 0, 0). Notice that this method has been used recently in a similar system in [START_REF] Neukirch | Survival of a recessive allele in a mendelian diploid model[END_REF]. In this paper the authors were interested in a system of diploid individuals, with dominance but without sexual preference and migration. Their center manifold was of dimension one.

Equilibria (ζ, 0, 0, 0, 0, ζ) and (0, 0, ζ, ζ, 0, 0): Those two equilibria with different monomorphic populations in the two patches are symmetrical. We only prove the results on (0, 0, ζ, ζ, 0, 0). The eigenvalues of the Jacobian Matrix at this fixed point are (0, -bβ

as well as the roots of the polynomial:

We were not able to determine all the polymorphic equilibria of the three dynamical systems governing the population dynamics in the two codominant and the dominant cases. Of course, thanks to our analysis of the cases without migration, we are able to give one unstable polymorphic equilibrium for each dynamical system:

• (δξ, ξ, δξ, δξ, ξ, δξ) for the first codominant model • (δ ξ , ξ , δ ξ , δ ξ , ξ , δ ξ ) for the second codominant model

, 1 for the model with dominance

To give an idea of the complexity of the dynamical systems under study, for the parameters b = 2, d = 1, c = 0.5, β = 1.1 and p = 5, Mathematica Inc. gives the numerical approximation of 9 polymorphic equilibria for the second codominant model.

Appendix E.2. A majority allele

Let us denote by ᾱ the complement of α ∈ {A, a}. In both codominant models, if z αα,i (0) ≥ z ᾱ ᾱ,i (0) for i = 1, 2, then z αα,i (t) ≥ z ᾱ ᾱ,i (t) for i = 1, 2 and for every positive t. Indeed if the (z αα ,i , αα ∈ {AA, Aa, aa}, i ∈ {1, 2}) evolve according to the dynamical system (A.1), we have

and if they evolve according to the dynamical system (A.2), we have

From these two equalities we find that the sets

are invariant under the dynamical systems (A.1) and (A.2).

Appendix F. The case of constant migration

For the sake of completeness, we study the case where all individuals have the same migration rate p, independent of the population state. The population dynamics is thus the solution to the following system of equations:

(F.1) We will show that in this case the equilibria are either null, or monomorphic, or have all their coordinates positive. More precisely Lemma 6. Let (z αα ,i , (α, α ) ∈ {A, a} 2 , i ∈ {1, 2}) be an equilibrium of Equation (F.1). Then, either z = 0, or all the coordinates of z are positive, or z has only two nonnull coordinates:

for an α ∈ {A, a}.

Let us denote by ᾱ the complement of α in {A, a}. The first step to get Lemma 6 is to prove the following Lemma.

Lemma 7. Let (z αα ,i , (α, α ) ∈ {A, a} 2 , i ∈ {1, 2}) be an equilibrium of Equation (F.1). Assume that there exists (α, i) ∈ {A, a} × {1, 2} such that z αα,i = 0. Then z = 0 or

and the other coordinates of z are null.

Proof. Let us first notice that 0 is an equilibrium and that if the population size of one patch is null, the population size of the other patch is also null. Thus we may assume in the rest of the proof without loss of generality that z AA,1 = 0, that N 1 > 0 and N 2 > 0. Then the first equation in (F.1) becomes:

and hence, z Aa,1 = z AA,2 = 0.

As z AA,2 = 0, using the same reasoning we get that z Aa,2 = 0. Now we are left with the last equation in (F.1), which can be written: 0 = (bβ -(d + cz aa,1 ))z aa,1 + p(z aa,2z aa,1 ). (F.2)

Similarly we get 0 = (bβ -(d + cz aa,2 ))z aa,2 + p(z aa,1z aa,2 ). (F.3) First, notice that a possible solution is z aa,1 = z aa 2 = (bβd)/c. Assume that z aa,1 z aa,2 . Then without loss of generality, we may assume that z aa,2 > z aa,1 . In this case, (F.2) entails: where R is the polynomial function which at X associates

This polynomial function is negative on I 1 := (0, (bβ-d)/c) and positive on I 2 := ((bβd)/c, ∞). Thus (F.5) implies that one of the z aa,i s' belongs to I 1 , and the other one to I 2 . This contradicts (F.4).

The second and last step to get Lemma 6 is to notice that if z is an equilibrium of (F.1) and if z Aa,1 = 0, then the second equation in (F.1) becomes: 2b N 1 z AA,1 z aa,1 + pz Aa,2 = 0.

This implies that z Aa,2 = 0, and z AA,1 = 0 or z aa,1 = 0. Then we may apply Lemma 7. This ends the proof of Lemma 6.