L Dekens
email: dekens@math.univ-lyon1.fr

S P Otto

V Calvez

The best of both worlds: combining population genetic and quantitative genetic models

Keywords: quantitative genetics, population genetics, heterogeneous environment, genetic architecture, polymorphism, asymptotic analysis

Numerous traits under migration-selection balance are shown to exhibit complex patterns of genetic architecture with large variance in effect sizes. However, the conditions under which such genetic architectures are stable have yet to be investigated, because studying the influence of a large number of small allelic effects on the maintenance of spatial polymorphism is mathematically challenging, due to the high complexity of the systems that arise. In particular, in the most simple case of a haploid population in a two-patch environment, while it is known from population genetics that polymorphism at a single major-effect locus is stable in the symmetric case, there exists no analytical predictions on how this polymorphism holds when a polygenic background also contributes to the trait. Here we propose to answer this question by introducing a new eco-evo methodology that allows us to take into account the combined contributions of a major-effect locus and of a quantitative background resulting from small-effect loci, where inheritance is encoded according to an extension to the infinitesimal model. In a regime of small variance contributed by the quantitative loci, we justify that traits are concentrated around the major alleles, according to a normal distribution, using new convex analysis arguments. This allows a reduction in the complexity of the system using a separation of time scales approach. We predict an undocumented phenomenon of loss of polymorphism at the major-effect locus despite strong selection for local adaptation, because the quantitative background slowly disrupts the rapidly established polymorphism at the major-effect locus, which is confirmed by individual-based simulations. Our study highlights how segregation of a quantitative background can greatly impact the dynamics of major-effect loci by provoking migrational meltdowns. We also provide a comprehensive toolbox designed to describe how to apply our method to more complex population genetic models.

Introduction

Biological motivation. Many species, if not most, evolve in heterogeneous habitats, where varying selection acts upon phenotypic traits in a manner that causes local adaptation. The genetic architecture that underlies those traits is known to present an array of possibilities, from major responses at one particular gene to diffuse polygenic reponses [START_REF] Slate | INVITED REVIEW: Quantitative trait locus mapping in natural populations: progress, caveats and future directions[END_REF]; Walsh and Lynch 2018). However, despite the boom in of genome sequencing of the last four decades, global conclusions on the conditions leading to a major gene or a polygenic response to local adaptation are yet to be drawn from empirical studies. For example, as reviewed in [START_REF] Walsh | Evolution and Selection of Quantitative Traits[END_REF], different conclusions on the genetic basis of the evolution of resistance to the insecticide BT toxin have emerged between field and lab experiments. Indeed, in the field, major-effects are more often found to be the main drivers of evolution of resistance, whereas a polygenic response is more commonly found in the lab (McKenzie and Batterham 1994), even if intensity of selection might not differ [START_REF] Groeters | Roles of Selection Intensity, Major Genes, and Minor Genes in Evolution of Insecticide Resistance[END_REF]. In more recent studies, divergent conclusions about the genetic basis of pathogen resistance in cattle have been reached in different regions of the world (major-effect in Australia: [START_REF] Turner | A genome-wide association study of tick burden and milk composition in cattle[END_REF], polygenic in the tropics: Porto-Neto et al. 2014). Other empirical studies also highlight cases where the genetic basis of local adaptation has a large variance in effect size, thus combining major and polygenic responses (see e.g. the stability of polymorphism at a major-effect locus underlying local adaptation in a symmetrical heterogeneous environment in the presence of a quantitative background due to a large number of small effect loci. Due to small perturbations induced by the quantitative component of the trait, polymorphism at the major-effect locus is lost both at low and high strengths of selection, below a certain level of migration. The first region of loss of polymorphism, at low selection intensities, is intuitively expected, as migration blends more strongly than selection differentiates. More surprising is the lost of polymorphism at high intensities of selection, where one would expect polymorphism at the major-effect locus to be strongly favoured. To our knowledge, this phenomenon, where quantitative differences displace polymorphism at a major-effect locus, has not yet been documented. We confirm that our analysis is consistent with individual-based simulations.

This case study suggests that the long-term influence of a quantitative polygenic background on the polymorphic equilibrium at major-effect loci can lead to unforeseen phenomena. In this work, we present an integrative framework that is meant to help analytically bridge population genetics and quantitative genetics. Our method goes deeper than previous models [START_REF] Lande | The response to selection on major and minor mutations affecting a metrical trait[END_REF]) by justifying in a certain regime of small variance that the traits are normally distributed around the major-effect alleles effects, thanks to new arguments of convex analysis. It allows a separation of time scales, which ultimately leads to the conditions for when the infinitesimal quantitative background slowly disrupts the rapidly established symmetrical polymorphism at the major-effect locus.

Furthermore, we provide a comprehensive toolbox that describes how to apply our methodology to more general cases in terms of number of major-effect loci, number of patches, and form of selection for haploid or diploid populations (see Appendices A and B).

Methods

Model

From a generic quantitative genetic model to a composite model.

We consider a haploid population reproducing sexually and characterized by a quantitative trait ζ in a heterogeneous environment with two habitats connected by constant migration at rate m 1 (from habitat 1 to habitat 2) and m 2 (from habitat 2 to habitat 1). Following classical models of quantitative genetics, we model each habitat i selecting toward a different optimum θ i with strength g i . Maladaptation and local uniform competition for resources (with intensity κ i in deme i) are sources of mortality leading to a per capita decline at rate:

-g i (ζ -θ i ) 2 -κ i N i ,
for individuals of trait ζ in habitat i (N i denotes the local population size). At time t ≥ 0, let n 1 (t, ζ) and n 2 (t, ζ) be the local trait densities in patches 1 and 2, and B[n i ](t, ζ) the number of individuals born with a trait ζ in habitat i, with reproduction occuring at rate λ i . The dynamics of the local trait densities read:

               ∂n1 ∂t (t, ζ) = λ 1 B[n 1 ](t, ζ) -g 1 (ζ -θ 1 ) 2 n 1 (t, ζ) -κ 1 N 1 (t) n 1 (t, ζ) +m 2 n 2 (t, ζ) -m 1 n 1 (t, ζ), ∂n2 ∂t (t, ζ) = λ 2 B[n 2 ](t, ζ) -g 2 (ζ -θ 2 ) 2 n 2 (t, ζ) -κ 2 N 2 (t) n 2 (t, ζ) +m 1 n 1 (t, ζ) -m 2 n 2 (t, ζ).
(1)

We can define the trait axis such that: θ := θ 2 = -θ 1 > 0 without loss of generality. We next describe the novel aspect of this work, which allows the trait ζ to be the sum of two components, a major-effect locus and a quantitative background z. We furthermore describe the sexual reproduction operator used.

major-effect. The first component comes from a locus where two alleles A/a are segregating. They have a major-effect on the trait: η A and η a . Inheritance of this locus is Mendelian.

Quantitative background. The second component, denoted by z ∈ R, represents the quantitative background due to infinitesimally small additive contributions to the trait from a large number of unlinked alleles. Although it comes from infinitesimally small contributions, z should not be thought of as being necessarily small, due to the large number of alleles contributing to it. We also assume that the major-effect locus is effectively unlinked with the small-effect ones.

Inheritance of the trait: an extension of the infinitesimal model. Let us recall that

the infinitesimal model, first introduced in [START_REF] Fisher | The Correlation between Relatives on the Supposition of Mendelian Inheritance[END_REF], provides a way to encode efficiently the inheritance of complex traits coming from a large number of alleles, each with small effects. The classical version states that an offspring receives a trait Z from its parents with traits Z 1 and Z 2 , where Z differs from the mean parental trait Z1+Z2 2 following a centered Gaussian law, with variance σ 2 2 . The latter accounts for the stochasticity of segregation, and therefore the variance is called the segregational variance. Specifically:

Z|Z 1 , Z 2 ∼ Z 1 + Z 2 2 + Y, Y ∼ N 0, σ 2 2 , Y ⊥ Z 1 , Z 2 .
The Mendelian view of the infinitesimal model has been discussed in [START_REF] Fisher | The Correlation between Relatives on the Supposition of Mendelian Inheritance[END_REF], [START_REF] Bulmer | The Effect of Selection on Genetic Variability[END_REF], and [START_REF] Lange | Central limit theorems of pedigrees[END_REF]: the common interpretation is that the trait results from a large number of small additive contributions at unlinked loci. For a more in depth description, see [START_REF] Barton | The infinitesimal model: Definition, derivation, and implications[END_REF].

Because the trait we are considering is a composite of a major-effect locus inherited according to Mendelian laws and an infinitesimal background, it is natural to use an extension of the infinitesimal model for this composite case. Now, the offspring's trait (A, Z) given their parents (A 1 , Z 1 ) and (A 2 , Z 2 ) reads:

(A, Z, ) | (A 1 , Z 1 ), (A 2 , Z 2 ) ∼ XA 1 + (1 -X)A 2 , Z 1 + Z 2 2 + Y , ( 2
)
where Y ∼ N 0, σ 2 2 follows a centered Gaussian law of variance σ 2 2 and X ∼ B 1 2 follows a Bernoulli law with parameter 1 2 (assuming fair meiosis). The random variables are independent of each other and of Z 1 , Z 2 , A 1 , A 2 .

Let us translate Eq. (2) into a continuous density model. Let n A i (z) (respectively n a i (z)) denote the density of individuals of patch i carrying allele A (respectively a) along with an infinitesimal background z, therefore having a trait ζ = η A + z (respectively, η a + z). In agreement with Eq. (2), the number of offspring born with the allele A and an infinitesimal contribution z in habitat i then reads:

B A σ [n A i , n a i ](z) = R 2 1 √ πσ exp - z -z1+z2 2 2 σ 2 × 1 N i n A i (z 1 ) n A i (z 2 ) + 1 2 n A i (z 1 ) n a i (z 2 ) + n a i (z 1 )n A i (z 2 ) dz 1 dz 2 = R 2 1 √ πσ exp - z -z1+z2 2 2 σ 2 n A i (z 1 ) n A i (z 2 ) + n a i (z 2 ) N i dz 1 dz 2 .
Similarly, the corresponding number of offspring born with the allele a and an infinitesimal part z reads:

B a σ [n A i , n a i ](z) = R 2 1 √ πσ exp - z -z1+z2 2 2 σ 2 × 1 N i n a i (z 1 ) n a i (z 2 ) + 1 2 n A i (z 1 ) n a i (z 2 ) + n a i (z 1 )n A i (z 2 ) dz 1 dz 2 = R 2 1 √ πσ exp - z -z1+z2 2 2 σ 2 n a i (z 1 ) n a i (z 2 ) + n A i (z 2 ) N i dz 1 dz 2 .
The operator reproduction B σ indicates that it is more relevant to model the dynamics of the two local allelic densities n a i , n A i , instead of n i (which is their sum). From now on, we will therefore adopt this point of view.

Remark 1: Bridging a population genetic model and a quantitative genetic model. Our model described above bridges the following population genetic and quantitative genetic models:

1. The one-locus haploid model in a two-patch environment, which considers two alleles A and a segregating at the same locus, each improving the survival chance in one of the habitats and being deleterious in the other. We recall that with symmetrical migration and selection, this model predicts that polymorphism at the focal locus is always stable, whenever the metapopulation persists (see Remark 4 and Proposition D.3 for a proof of this fact).

2. The quantitative genetic model from [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF], which studies the eco-evo dynamics of a quantitative trait in a heterogeneous environment, where the trait is inherited according to the standard version of the infinitesimal model. Our work can be seen as an extension of this model, to which we add the segregation of two major-effect alleles at a single locus. Moreover, one can notice that if one major-effect allele fixes (loss of polymorphism), the two models are equivalent. Because [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF] gives a complete analytical description of the outcomes of their system (in the small segregation variance regime), the outcomes for our present study are known given the fixation of a major-effect allele. Therefore, our study focuses on the description of polymorphism at the major-effect locus and its stability.

Dimensionless system.

Let us rescale Eq. (1) according to:

η A := η A θ , z := z θ , g i := g i θ 2 λ 1 , m i := m i λ 1 , ε := σ θ , t := ε 2 λ 1 t, , α := κ 1 κ 2 , λ := λ 2 λ 1 ,
and introduce the rescaled trait densities:

n A ε,i (t, z) := κ i λ 1 n A i (t, z), n a ε,i (t, z) := κ i λ 1 n a i (t, z).
so that Eq. (1) reads:

                               ε 2 ∂n A ε,1 ∂t (t, z) = B ε A (n A ε,1 , n a ε,1 )(t, z) -g 1 (z + η A + 1) 2 n A ε,1 (t, z) -N ε,1 (t) n A ε,1 (t, z) +α m 2 n A ε,2 (t, z) -m 1 n A ε1 (t, z), ε 2 ∂n a ε,1 ∂t (t, z) = B ε a (n a ε,1 , n A ε,1 )(t, z) -g 1 (z + η a + 1) 2 n a ε,1 (t, z) -N ε,1 (t) n a ε,1 (t, z) +α m 2 n a ε,2 (t, z) -m 1 n a ε1 (t, z), ε 2 ∂n A ε,2 ∂t (t, z) = λ B ε A (n A ε,2 , n a ε,2 )(t, z) -g 2 (z + η A -1) 2 n A ε,2 (t, z) -N ε,2 (t) n A ε,2 (t, z) + m1 α n A ε,1 (t, z) -m 2 n A ε2 (t, z), ε 2 ∂n a ε,2 ∂t (t, z) = λ B ε a (n a ε,2 , n A ε,2 )(t, z) -g 2 (z + η a -1) 2 n a ε,2 (t, z) -N ε,2 (t) n a ε,2 (t, z) + m1 α n a ε,1 (t, z) -m 2 n a ε2 (t, z), (3)
where the rescaled reproduction operator is given by:

B A ε (n A ε,i , n a ε,i )(t, z) = 1 √ πε R 2 exp -(z -z1+z2 2 ) 2 ε 2 n A ε,i (t, z 1 ) n A ε,i (t, z 2 ) + n a ε,i (t, z 2 ) N ε,i (t) dz 1 dz 2 . (4)
2.2 Derivation of a moment-based system in the regime of small variance ε 2 ≪ 1

In this subsection, we explain how we derive a closed moment-based ODE system on which the separation of time scale analysis will be conducted, starting from the PDE system (3) based on the trait distributions, in the regime of small variance ε 2 ≪ 1. To do so, we justify that the quantitative background values are approximately normally distributed among bearers of the same major-effect allele. Moreover, the mean of these quantitative background values is the same for individuals in the same patch. This implies in particular that the main driver for trait divergence within each habitat is the major-effect locus.

Gaussian approximations of quantitative background values in the regime of small variance: a formal analysis.

We choose to place our study in a regime where the amount of diversity introduced by the segregation of the infinitesimal background at each event of reproduction is small in comparison to the difference between the habitats' optima:

σ 2 θ 2 ≪ 1 =⇒ ε 2 ≪ 1.
In this regime of small variance, the trait distributions are expected to converge to Dirac masses. Our focus is to give an accurate description of the distribution near this limit. To do so, we extend a small variance methodology introduced by Diekmann As the analytical crux heavily relates to the singular nature of the trait distributions n ε as Dirac masses, the method consists in defining proxies U ε from n ε through a suitable transformation so that such proxies are regular functions (by comparison to Dirac masses) and their asymptotic analysis is easier. Studying them often induces a reduction in the complexity of the system while still retaining fundamental quantitative information about the distributions, such as around which traits they are concentrated. Here, we follow quantitative genetic studies that use the infinitesimal model according to the same methodology (Garnier et

U A ε,i (resp. U a ε,i ): n A ε,i = 1 √ 2πε e - U A ε,i ε 2 , n a ε,i = 1 √ 2πε e - U a ε,i ε 2 . ( 5
)
A helpful analogy is to take the example of a spiky Gaussian distribution with small variance ε 2 for n A ε,i . Then U A ε,i is a smooth quadratic function (even when ε ≪ 1). Figure 1 displays an example of this kind of exponential transformation (called Hopf-Cole transformation in scalar conservation laws).

A key observation to deduce the traits around which the distribution concentrates is that it does so at the minima (zero) of U ε . As the proxies U A ε and U a ε are expected to be more regular in the regime of small variance, they are thought to be the right object on which to perform a Taylor expansion series to gain information on the asymptotic distributions in the limit of small variance (see [START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF]. We therefore define u A 0,i (resp. u a 0,i ) as the leading term in the Taylor expansion of

U A ε,i (resp. U a ε,i ) : U A ε,i = u A 0,i + ε 2 u A 1,i + ε 4 v A ε,i , U a ε,i = u a 0,i + ε 2 u a 1,i + ε 4 v a ε,i (6)
where u A 1,i and u a 1,i are the next term in the Taylor expansion, and ε 4 v A ε,i and ε 4 v a ε,i are the residues. Calvez, Garnier, and Patout (2019) provides the tools to control these residues and thus rigorously justify that (6) is an admissible Taylor expansion; adapting them is left for future work.

Characterization of the main terms u A

0,i and u a 0,i to justify Gaussian approximations.

The first step of the analysis in the regime of small variance is the characterization of the main terms u A 0,i and u a 0,i . Indeed, in the regime of small variance, these have to satisfy a strong constraint that arises naturally for the contribution of the infinitesimal model reproduction operator term to remain wellbalanced within (3). In [START_REF] Garnier | Adaptation to a changing environment: what me Normal?[END_REF] and [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF], where the standard infinitesimal model operator is used, this constraint yields the analogous main term to be quadratic, which implies that the trait distribution is approximately Gaussian, with a small variance ε 2 . However, here, the arguments given in [START_REF] Garnier | Adaptation to a changing environment: what me Normal?[END_REF] and used in [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF] are not sufficient, due to the mixing of alleles between patches and the discrete nature of Mendelian inheritance. However, we extend the convex analysis to circumvent this limitation (Proposition 2.1) and identify u A 0,i and u a 0,i as the same quadratic function z →

(z-z * i ) 2 2
, where z * i ∈ R is to be determined later in the analysis. Assuming that ( 6) is an admissible Taylor expansion (which is suggested by the analysis of Calvez, Garnier, and Patout (2019)), this result is crucial as it justifies the following formal Gaussian approximations of n A ε,i and n a ε,i (i ∈ {1, 2}):

n A ε,i (z) = e - -(z-z * i ) 2 2ε 2 √ 2πε e -u A 1,i (z)+O(ε 2 ) , n a ε,i (z) = e - -(z-z * i ) 2 2ε 2 √ 2πε e -u a 1,i (z)+O(ε 2 ) . ( 7
)
Hence, to the leading order, n A ε,i and n a ε,i are formally Gaussian, centered at the same quantitative contribution z * i , with the same variance ε 2 . However, they differ in the next-order, which involves the corrector terms u A 1,i and u a 1,i , which generate asymmetries in the distributions. To support (7), we first derive the following constraints (C) on the main terms u A 0,i and u a 0,i . In order for the contribution of both reproduction operators B A ε and B a ε to remain well-balanced with the other biological phenomena in the regime of small variance in (3), u A 0,i and u a 0,i formally need to satisfy the following (see Appendix C for the details):

                   ∀z ∈ R, max sup z1,z2 u A 0,i (z) -z -z1+z2 2 2 -u A 0,i (z 1 ) -u A 0,i (z 2 ), sup z1,z2 u A 0,i (z) -z -z1+z2 2 2 -u A 0,i (z 1 ) -u a 0,i (z 2 ) = 0, ∀z ∈ R, max sup z1,z2 u a 0,i (z) -z -z1+z2 2 2 -u a 0,i (z 1 ) -u a 0,i (z 2 ), sup z1,z2 u a 0,i (z) -z -z1+z2 2 2 -u A 0,i (z 1 ) -u a 0,i (z 2 ) = 0. (C)
We next state the following proposition, which characterizes the main terms u A 0,i and u a 0,i as aforementioned.

Proposition 2.1. Let u A 0 and u a 0 satisfying Eq. (C) positive almost everywhere and cancelling somewhere. Then, there exists z * ∈ R such that:

∀z ∈ R, u A 0 (z) = u a 0 (z) = (z -z * ) 2 2 . ( 8
)
The conditions on u A 0,i and u a 0,i in Proposition 2.1 (positive everywhere and cancelling somewhere) are explained in Appendix C. In the Appendix B, we actually state and prove a stronger result Proposition B.1, which generalizes Proposition 2.1 to more complex genetic architectures.

Consequently, assuming that (6) is the correct ansatz so that we can control the residues in (7) (which the analysis of Calvez, Garnier, and [START_REF] Calvez | Asymptotic analysis of a quantitative genetics model with nonlinear integral operator[END_REF] suggests and provides a framework to show), using the result of Proposition 2.1 in (5) and (6) leads to (7).

Moment-based system in the regime of small variance

This section follows directly the results of the previous one, where we showed formally that, in each habitat, the two allelic trait distributions n A ε,i and n a ε,i can be approximated by the same Gaussian distribution. We present here how the latter allows us to close the moment-based system obtained from integrating (3).

First, we derive formal expansions of the first moments (population size, mean trait, variance and skew) of n A ε,i and n a ε,i when ε 2 ≪ 1, thanks to (6) and (7) (as in Dekens 2022):

                       N A ε,i := R n A ε,i (z) dz = e -u A 1,i (z * i ) 1 + ε 2 ∂z u A 1,i (z * i ) 2 2 - ∂zz u A 1,i (z * i ) 2 -v A i,ε (z * i ) + O(ε 4 ), z A ε,i := R z n A ε,i (z) N A ε,i dz = z * i -ε 2 ∂zu A 1,i (z * i ) + O(ε 4 ), σ A ε,i 2 := R (z -z A ε,i ) 2 n A ε,i (z) N A ε,i dz = ε 2 + O(ε 4 ), ψ A ε,i 3 := R (z -z A ε,i ) 3 n A ε,i (z) N A ε,i dz = O(ε 4 ). ( 9
)
Using (9) when integrating (3), we can close the infinite system of moments in the regime of small variance, producing a system of eight ODEs governing the dynamics of the four allelic subpopulation sizes N a ε,1 , N A ε,1 , N a ε,2 , N A ε,2 and the four allelic local mean quantitative traits

z a ε,1 , z A ε,1 , z a ε,2 , z A ε,2 :                                                                                                    ε 2 d N a ε,1 dt = N a ε,1 -N A ε,1 + N a ε,1 N a ε,1 -g 1 z a ε,1 + η a + 1 2 N a ε,1 + α m 2 N a ε,2 -m 1 N a ε,1 , +O(ε 2 ), ε 2 d N A ε,1 dt = N A ε,1 -N A ε,1 + N a ε,1 N A ε,1 -g 1 z A ε,1 + η A + 1 2 N A ε,1 + α m 2 N A ε,2 -m 1 N A ε,1 +O(ε 2 ), ε 2 d N a ε,2 dt = λ N a ε,2 -N A ε,2 + N a ε,2 N a ε,2 -g 2 z a ε,2 + η a -1 2 N a ε,2 + m1 α N a ε,1 -m 2 N a ε,2 +O(ε 2 ), ε 2 d N A ε,2 dt = λ N A ε,2 -N A ε,2 + N a ε,2 N A ε,2 -g 2 z A ε,2 + η A -1 2 N A ε,2 + m1 α N a ε,1 -m 2 N a ε,2 +O(ε 2 ), ε 2 d z a ε,1 dt = ε 2 2g 1 -1 -η a -z a ε,1 + z A ε,1 -z a ε,1 2 N A ε,1 Nε,1 + α m 2 N a ε,2 N a ε,1 z a ε,2 -z a ε,1 + O(ε 4 ), ε 2 d z A ε,1 dt = ε 2 2g 1 -1 -η A -z A ε,1 + z a ε,1 -z A ε,1 2 N a ε,1 Nε,1 + α m 2 N A ε,2 N A ε,1 z A ε,2 -z A ε,1 + O(ε 4 ), ε 2 d z a ε,2 dt = ε 2 2g 2 1 -η a -z a ε,2 + z A ε,2 -z a ε,2 2 N A ε,2 Nε,2 + m1 α N a ε,1 N a ε,2 z a ε,1 -z a ε,2 + O(ε 4 ), ε 2 d z A ε,2 dt = ε 2 2g 2 1 -η A -z A ε,2 + z a ε,2 -z A ε,2 2 N a ε,2 Nε,2 + m1 α N A ε,1 N A ε,2 z A ε,1 -z A ε,2 + O(ε 4 ). ( 10
)
Biological description of the equations of the moment-based system (10). The first four equations encoding the dynamics of the allelic subpopulations sizes involve four terms, that we describe using the first equation for N a ε,1 . The first term N a ε,1 is a growth term, the second one -(N A ε,1 + N a ε,1 )N a ε,1 is a non-linear negative death term by competition, proportional to the total subpopulation size. The third one -g 1 z a ε,1 + η a + 1 2 N a ε,1 is a negative death term by selection (with strength g 1 ),

which is more lethal when the allelic local mean trait z a ε,1 + η a is far from the local optimum -1. The last migration term α m 2 N a ε,2 -m 1 N a ε,1 represents the asymmetrical transfer of populations between the two patches.

The last four equations encoding the dynamics of the allelic local mean quantitative traits involve three different terms that we describe, taking for reference the first of these equations for z a ε,1 . The first term is the selection gradient that pushes the total mean trait z a ε,1 + η a towards the local optimum -1, with an intensity 2ε 2 g 1 , proportional to the intensity of selection g i and the small variance of the quantitative trait ε 2 (in agreement with the Gaussian approximation (7)). The second term

z A ε,1 -z a ε,1 2 N A ε,1
Nε,1 does not exist in the analogous moment-based system in Dekens (2022) (without the major-effect locus), as it originates from the segregation of A/a at the major-effect locus. It describes a force which pushes each allelic mean quantitative component towards one another within the same habitat due to the mixing effect of the infinitesimal model. It is consistent with the result provided by Proposition 2.1 and the Gaussian approximations (7), which are centered at the same quantitative component z * i , close to both z A ε,i and z a ε,i according to the second line of the expansions (9). The last term α m 2

N a ε,2 N a ε,1 z a ε,2 -z a ε,1
relates to the effect of the transfer of population by migration onto the mean quantitative component: it pushes the local mean quantitative components corresponding to the same major-effect allele z a ε,2 and z a ε,1 towards one another. Remark 2: Selection shifts the allelic local mean quantitative trait slowly. In the last four equations of (10), there is a noticeable difference between the first term, proportional to ε 2 , and the other two terms, which are of order 1. This demonstrates the fact that, in the regime of small variance, selection shifts the local mean quantitative traits very slowly toward the local optima compared to how fast the other two terms intervene in the equation (describing selection on the major-effects alleles and migration). Notice also that the time scale in which the differential system (10) is written (ε 2 d• dt ) is the correct one to capture this slow shift.

Remark 3: Magnitude of the residues in (10). In the system (10), the difference in the system between the residues in the first four equations on the local sizes of population of order O(ε 2 ) and the ones in the last four equations on the mean quantitative components of order O(ε 4 ) is consistent with the analysis of Patout (2020) (see in particular Theorem 1.4).

Separation of time scales: slow-fast analysis

As highlighted by Remark 2, the shift of allelic local mean quantitative components z A ε,i and z a ε,i occurs on a slower time scale than growth, death and transfer of populations for the allelic subpopulation sizes (first four equations of (10)) and than the two relaxing forces of gene flow and segregation for the allelic local mean quantitative traits (last two terms of the last four equations of (10)). Therefore, in this subsection, we show that the moment-based system (10) has a particular structure (up to a change in variables) that allows the possibility to separate two different time scales, which can be interpreted as fast ecological time scale (including selection on the major-effects locus) and slow quantitative evolutionary time scales.

First, we need to transform (10) into an equivalent system which has a suitable form to prove the separation of time scales. This requires the following change of variables, which is motivated by the formal analysis of Section 2.2 (especially the results of Proposition 2.1):

δ a ε = z a ε,2 -z a ε,1 2ε 2 , δ A ε = z A ε,2 -z A ε,1 2ε 2 , δ ε = z A ε,1 + z A ε,2 -z a ε,1 -z a ε,2 4ε 2 , Z ε = z A ε,1 + z A ε,2 + z a ε,1 + z a ε,2 4 . ( 11
) Z ε can be interpreted as the mean infinitesimal part of the metapopulation, δ ε the spatial average of the local difference between the two allelic mean infinitesimal parts, δ A ε and δ a ε the equivalent term among bearers of A and a, respectively (see an illustration of those new variables in Fig. 2). The quantities defining δ ε , δ A ε , δ a ε are divided by ε 2 because Remark 2 suggests that z A ε,1 , z a ε,1 , z A ε,2 and z a ε,2 all relax quickly towards the same value due to the fast action of gene flow and segregation, with an error of order ε 2 . δ ε ,δ A ε and δ a ε (in red), introduced in (11). This figure displays a situation where the two major-effect alleles are segregating in both habitats in a symmetrical fashion. The graph represents the two local trait densities for each of the two alleles: ñA 1,ε , ñA 2,ε , ña 1,ε , ña 2,ε (the same color is for the same deme, and the same linestyle is for the same major-effect allele), as a function of the trait ζ = z + η A (resp. z + η a ), where z is the infinitesimal contribution and η A (resp η a ) is the effect of the major-effect allele. In red, we indicate a visualization of the new variables introduced in (11). Z ε is the mean infinitesimal part of the metapopulation, δ ε the spatial average of the local difference between the two allelic mean infinitesimal parts, δ A ε and δ a ε the spatial discrepancies in the mean infinitesimal parts per allele. Note the difference in notation between the trait densities ñA i,ε and the infinitesimal contribution densities n A i,ε (which are the ones used in the analysis), which are linked by

ε 2 δ a ε ε 2 δ A ε ε 2 δ ε + 2η Z ε
n A i,ε (z) = ñA i,ε (z + η A ) (respectively z + η a for ña i,ε ).
After the change in variables (11), we obtain a new system (shown in its explicit form Eq. (28) in Appendix D) which can be written compactly as follows

ε 2 d Ȳε dt = G( Ȳε , Z ε ) + ε 2 ν G ε (t), dZε dt = -(g 1 + g 2 ) Z ε + F ( Ȳε ) + ε 2 ν F ε (t), (P ε )
where

Ȳε := N a 1,ε , N a 1,ε , N A 1,ε , N A 2,ε , δ a ε , δ A ε , δ ε denotes the vector of fast variables, located in a set denoted Ω := R * + 4 × R 3 . The two smooth functions G ∈ C ∞ (Ω × R) and F ∈ C ∞ (Ω) encode
respectively the fast and slow dynamics. Moreover, the functions ν G ε and ν F ε are residues that are uniformly bounded w.r.t ε.

Biological interpretation of the different timescales in (P ε ). As the first four coordinates

of the fast variable Ȳε are the allelic subpopulation sizes, the function G(•, Z ε ) describes the fast dynamics of growth, death and transfer of populations occurring when the mean quantitative component is at the value Z ε . The fast timescale of the dynamics of Ȳε can be interpreted as the ecological time scale (including selection on the major-effects locus). On the contrary, the dynamics of the slow variable Z ε , which is the mean quantitative component, are driven by the shift by selection -(g 1 +g 2 )Z ε and the demographic feedback F ( Ȳε ), on a slower timescale, which we interpret as the quantitative evolutionary time scale. Indeed, notice that the time derivatives are different between the two lines of P ε : the first line involves ε 2 d• dt , whereas the small factor ε 2 is absent in the second line.

Convergence to a simplified limit system. The slow-fast analysis developed in Appendix D is dedicated to show that, when ε goes to 0, the solutions of (P ε ) converge to the solutions of the following limit system which separates the ecological and evolutionary time scales

G( Ȳ , Z) = 0, dZ dt = -(g 1 + g 2 ) Z + F ( Ȳ ). (P 0 )
The first line of (P 0 ) is an algebraic system defining the slow manifold, constituted by the fast ecological equilibria Ȳ corresponding to a value Z of the evolutionary variable (these are formally defined by

{ Ȳ ∈ Ω, such that G( Ȳ , Z) = 0}
). The second line describes the dynamics of the slow variable Z constrained to occur on the slow manifold. The convergence result linking (P ε ) to (P 0 ) is stated by the following:

Theorem 2.1. For ( Ȳ , Z) a solution of (P 0 ), there exists T * > 0 such that, for 0 < ε < 1, any solution ( Ȳε , Z ε ) of (P ε ) on [0, T * ] converges to ( Ȳ , Z) uniformly on [0, T * ], as ε goes to 0 and ( Ȳε (0), Z ε (0)) goes to ( Ȳ (0), Z(0)).

The proof the Theorem 2.1 is detailed in Appendix D. The main argument relies crucially on the stability of the fast equilibria at any level defined by a value of the slow variable Z ∈] -1, 1[ (Proposition D.3, Proposition D.4), ensuring that, at the limit, the fast dynamics converge quickly toward the slow manifold and not away from it. The stability argument is completed by the algebraic description of the slow manifold: we show that, for every level Z ∈]-1, 1[, there exists a single ecological equilibria Ȳ satisfying G( Ȳ , Z) = 0 (Proposition D.1, Proposition D.2). We also summarize in Fig. 7 the links between the different systems, propositions and theorem involved in the slow-fast analysis.

Remark 4: The one-locus haploid model's equilibrium is part of the fast equilibrium corresponding to the level Z = 0. The one-locus haploid model is equivalent to the first four differential equations of (10) on the allelic sizes of each subpopulation, with (z a ε,1 , z a ε,1 , z a ε,1 , z a ε,1 ) = (0, 0, 0, 0) (no infinitesimal part -we can obtain from these equations a system describing the allelic frequencies and local population sizes

(p 1 , p 2 , N 1 , N 2 ) := N A 1 N A 1 +N a 1 , N A 2 N A 2 +N a 2 , N A 1 + N a 1 , N A 2 + N a 2 ,

dropping the ε that is a parameter of the infinitesimal part). Applying Proposition D.1 with Z = 0 gives a unique equilibrium satisfying the first four equations, which is the one found with the one-locus haploid model. One can thus interpret the symmetrical polymorphic equilibrium of the one-locus haploid model as a fast equilibrium in the model presented in this article. It is therefore stable (Proposition D.3) whenever it entails positive population sizes (same condition as in Proposition D.1).

Remark 5: Degrees of freedom of the slow manifold compared to Dekens 2022. Proposition D.1 states that for every level Z ∈] -1, 1[, there exists a single fast equilibrium Ȳ such that G( Ȳ , Z) = 0. This implies that there are fewer degrees of freedom in the subsystem (S 0 (Z)) defining the four allelic subpopulations sizes (see the details in Appendix D) than in the analogous system of two equations from the analysis done in [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF], that can be obtained in the case where one allele has fixed (up to a translation). Indeed, [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF] shows that the analogous system can have up to three algebraic solutions depending on the parameters. The result of Proposition D.1 is thus unexpected, since S 0 (Z) has twice the number of equations and variables.

3 Results: stability of polymorphism at the major-effect locus in the limit system

This section follows naturally the separation of timescales shown in Section 2.3 and focuses on the study of the stability of polymorphism at the major-effect locus in the limit system (P 0 ), in the presence of a quantitative background contributing additively to the trait under selection. To be able to derive analytical conditions, we assume henceforth a symmetrical environment setting (in migration rates, selection strengths, carrying capacities, reproduction rates and major-effect allelic effects):

m := m 1 = m 2 , g := g 1 = g 2 , α = 1, λ = 1, η := η A = -η a > 0.
Under these symmetrical conditions, and in the absence of any quantitative background, we recall that there exists a symmetrical polymorphic equilibrium in the one-locus haploid model, which is always stable (see Proposition D.3 for a proof). This symmetrical polymorphic equilibrium in the one-locus model corresponds, in our model which considers additionally the additive contribution of a quantitative background on the trait, to the fast equilibrium Y * associated to the level Z * = 0 (Z * = 0 corresponds to the average quantitative trait between patches cancelling). Because the property of the fast equilibrium does not necessarily transpose to a global equilibrium over multiple timescales, we are therefore interested in the following questions:

1. Does a symmetrical polymorphic equilibrium for the global limit system (P 0 ) exist at the level Z * = 0, ie: does the pair of variables (Z * , Ȳ * ) defined above cancel both the first line and the right-hand side of the second line of (P 0 )?

2. When the symmetrical polymorphic equilibrium exists, is it always stable ? Or, in the long-run, can the slowly evolving infinitesimal background undermine the rapidly established polymorphism at the major-effect locus, even though the latter is appears favored for local adaptation?

3. If so, can our analysis predict in which range of parameters of migration rate m, selection strength g and major-effect η does that phenomenon occur?

In a first part, we present the results of our analysis to answer these questions. We also provide illustrations of the complex patterns that can emerge in terms of parameters range, as the studied phenomenon of disturbance of the polymorphism at the major-effect locus by the infinitesimal background exhibits non-monotonic behaviours with regard to each parameter.

In a second part, we confirm the results of the first part thanks to individual-based simulations.

Analytical predictions

The results of this section indicate that the unconditional stability of the polymorphism in the OLM can be disturbed by the presence of a quantitative background, for a substantial range of parameters, including, surprisingly, at the strongest selection levels. The interpretation of Remark 4 offers the idea that the infinitesimal background slowly disrupts the rapidly established symmetrical polymorphism at the major-effect locus.

Existence of a symmetrical polymorphic equilibrium. We first show that a symmetrical polymorphic equilibrium can exist under a range of parameters specified in Proposition 3.1, as a stationary state of the limit system (P 0 ), hence a solution of the explicit version of the latter:

                                             N a 1 -N A 1 + N a 1 N a 1 -g [Z -η + 1] 2 N a 1 + m(N a 2 -N a 1 ) = 0, N A 1 -N A 1 + N a 1 N A 1 -g [Z + η + 1] 2 N A i + m(N A 2 -N A 1 ) = 0, N a 2 -N A 2 + N a 2 N a 2 -g [Z -η -1] 2 N a 2 + m(N a 1 -N a 2 ) = 0, N A 2 -N A 2 + N a 2 N A 2 -g [Z + η -1] 2 N A 2 + m(N A 1 -N A 2 ) = 0, 2g -m δ a N a 2 N a 1 + N a 1 N a 2 + δ A -δ a 4 N A 2 N a 2 +N A 2 + N A 1 N a 1 +N A 1 + δ 2 N A 2 N a 2 +N A 2 - N A 1 N a 1 +N A 1 = 0, 2g -m δ A N A 2 N A 1 + N A 1 N A 2 + δ a -δ A 4 N a 2 N a 2 +N A 2 + N a 1 N a 1 +N A 1 + δ 2 N A 1 N a 2 +N A 2 - N a 2 N a 1 +N A 1 = 0, -δ 2 -2 g η + m δ A 2 N A 2 N A 1 - N A 1 N A 2 -δ a 2 N a 2 N a 1 - N a 1 N a 2 = 0, -2 g Z + m δ a 2 N a 2 N a 1 - N a 1 N a 2 + δ A 2 N A 2 N A 1 - N A 1 N A 2 + δ A -δ a 4 N A 2 N A 2 +N a 2 - N A 1 N A 1 +N a 1 + δ 2 N A 1 N A 1 +N a 1 + N A 2 N A 2 +N a 2 -1 = 0. ( 12
)
Proposition 3.1. There exists a unique polymorphic equilibrium corresponding to the infinitesimal average Z = 0 under the condition:

g(η 2 + 1) < 1 ∨ m < 2 g 2 η 2 g(η 2 + 1) -1 -g(η 2 + 1) + 1 . ( 13
)
The allelic local population sizes corresponding to this equilibrium satisfy the property:

N a, * 1 = N A, * 2 , N a, * 2 = N A, * 1 ,
and for both alleles, the spatial discrepancies between the mean infinitesimal parts of the two patches per allele are the same:

δ A, * = δ a, * .
Therefore, this polymorphic equilibrium at Z * = 0 is called symmetrical.

The proof uses the results of Proposition D.1 and Proposition D.2 and is shown in Appendix E.

Stability of the symmetrical polymorphic equilibrium. Let us recall the limit system

(P 0 ): G( Ȳ , Z) = 0, dZ dt = -2 g Z + F ( Ȳ ). ( 14
)
The stability of the symmetrical polymorphic equilibrium denoted (0, Ȳ * ) described above in Proposition 3.1 is studied in the same manner as in [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF]. Because the differential equation on Z in the second line of (P 0 ) involves both Z and Ȳ , it is not sufficient to do a standard linear analysis. The first step is to express the solution to G(•, Z) = 0 as a function of Z: Ȳ (Z). This is possible thanks to the implicit function theorem used in the vicinity of the symmetrical polymorphic equilibrium, because

[∂ Ȳ G] | Z=0, Ȳ = Ȳ * is invertible (thanks toProposition D.
3 and Proposition D.4). This step allows us to recast (P 0 ) as

dZ dt = F(Z) := -2 g Z + F Ȳ (Z) , Z ∈] -1, 1[. ( 15
)
The stability of the symmetrical polymorphic equilibrium (0, Ȳ * ) can now be analysed by using the chain rule of differentiation on the right-hand side of (15). We obtain that the symmetrical polymorphic equilibrium is asymptotically locally stable if and only if

0 < 2g + ∂ Ȳ F • [∂ Ȳ G] -1 ∂ Z G Z=0, Ȳ = Ȳ * .
Due to the large number of dimensions involved, the explicit formula of the latter is too long to be given here. The patterns resulting from the numerical analysis of the stability of the symmetrical polymorphic equilibrium for four values of the effect of the major-effect locus η ∈ {0.5, 0.7, 1, 1.3} are computed in Fig. 3. For each value of η, the region of the stability of the polymorphism is indicated in yellow with selection strength (g, x-axis) and migration rate (m, y-axis) varying in [0, 3]. We can first observe that these yellow regions have complex boundaries, and exhibit non-monotonic behaviours with regard to both migration rate m and selection strength g. These are not predicted by the one-locus symmetric model (OLM), which states that polymorphism is maintained everywhere under the dashed yellow line, which represents the extinction threshold without the quantitative component (computed thanks to the viability condition (13) stated in Proposition D.1 for Z = 0). The latter leads to the conclusion that, when it occurs, the instability of polymorphism at the major-effect locus shown by our analysis results stems from the presence of the quantitative background due to small-effect loci.

To describe the non-monotonic behaviour with respect to increasing selection strengths, one can consider holding a constant intermediate migration rate and increase selection (going left to right on a horizontal line in Fig. 3a, Fig. 3b, Fig. 3c and Fig. 3d). While the polymorphism at the major-effect locus is not stable with weak selection, stability is gained at an intermediate level of selection that depends on the migration rate and subsequently lost at a higher level of selection. This non-monotonic behaviour when increasing selection levels is quite robust with regard to different values of η, as shown by the different panels in Fig. 3 (even if the effect is attenuated when η = 1 in Fig. 3c, which means that the major-effects coincide with the local optima). When selection is weak compared to migration (left sides of Fig. 3a, Fig. 3b, Fig. 3c and Fig. 3d), it is expected that the relative blending by migration, which is strong compared to the divergent forces of local selection, provokes the loss of polymorphism. The loss of polymorphism at the major-effect locus is more surprising and counter-intuitive as one would expect that the bonus provided by polymorphism at the major-effect locus, which helps subpopulations to be locally adapted, would be even more important at stronger selection levels, and therefore maintained. Unfortunately, the explicit mathematical expression of ( 15) is too involved to be truly informative about what is the cause of the loss of polymorphism at the major-effect locus with strong selection. We recommend the reader interested in this to consult the next section presenting the results of individual-based simulations, which provide insights on the origin of this phenomenon.

Pushing further the numerical analysis of polymorphic equilibria. Here, we show a numerical analysis of all the equilibria of the limit system (P 0 ) in Fig. 4. To do so, we use the autonomous differential equation (15) derived previously thanks to the implicit function theorem (used on the whole interval Z ∈]-1, 1[ thanks to Proposition D.3 and Proposition D.4). From (15), (Z, Ȳ (Z)) is a polymorphic equilibrium if F(Z) = 0, and this equilibrium is locally stable if F ′ (Z) < 0.

Even if the complexity of the limit system is still too great to be analytically solved (due to the implicit nature of the function Ȳ defined by the relation G( Ȳ (Z), Z) = 0), we show in Fig. 4 the phase lines corresponding to the limit equation ( 15), when the migration rate and the effect size of the major-effect locus are held constant (m = 0.1, η = 0.5) and the selection strength varies (the lighter the color, the stronger the selection). Solid lines indicate that the system is polymorphic, whereas dotted lines indicate that one major-effect allele has fixed. Every intersection of the zero horizontal line and a solid colored line with a negative slope indicates a locally stable polymorphic equilibrium (conversely, a positive slope indicates an unstable equilibrium).

This figure is consistent with the analysis of Section 3.1 and Fig. 3a: at Z = 0, all the curves return to 0 (F (Z)), confirming that a polymorphic equilibrium exists when the mean contribution of the smalleffect loci is 0. Their local slope indicates the stability of this equilibrium (stable if negative, unstable if positive). Furthermore, Fig. 4 gives insights on the existence of asymmetrical polymorphic equilibria. Particularly, it seems that such equilibria exist for a narrow window of intermediate selection strength: the green curve corresponding to g = 0.86 displays two mirrored stable asymmetrical polymorphic equilibria at Z ≈ ±0.5 (indicated by the red arrows), which is hard to predict analytically due to the high orders of polynomials involved. Moreover, such equilibria are presumably quite subtle to catch in individual-based simulations, because the window of selection and the basin of attraction are both narrow. However, this illustrates the new and unsuspected insights that can be obtained from this composite model. ) and the selection strength varies (the lighter the color, the stronger the selection). Solid curves indicate that the system is polymorphic, whereas dashed curves indicate that one major-effect allele has fixed. Every intersection of the horizontal black line and a solid colored curve with a negative (resp. positive) slope indicates a locally stable (resp. unstable) polymorphic equilibrium. The darker curve with weak selection g = 0.01 has a positive slope at Z = 0 (unstable), the following curves have a negative slope at Z = 0 (stable for selection between g = 0.18 and g = 0.86), and finally the lightest curves have a positive slope at Z = 0 (unstable for g ≥ 1.03), which is consistent with Fig. 3a. Note that there exists additionally two mirrored asymmetrical polymorphic equilibria for g = 0.86, for Z ≈ ±0.5 (indicated by the red arrows), which were unsuspected prior to this numerical analysis.

Individual-based simulations

In this part, we confirm the results given by our analysis on the stability of the symmetrical polymorphic equilibrium, using individual-based simulations conducted with the software SLiM (Haller and Messer 2019). We focus on the gain and loss of polymorphism with regard to increasing selection, when η ̸ = 1, for symmetrical and asymmetrical initial conditions. For each set of parameters, we ran 20 replicate simulations. The results for the major-effect locus are displayed in Fig. 5 and Fig. 9, both with a quantitative background (left panel) and without (right panel). The simulations confirm that variation is maintained only for intermediate levels of selection (as measured by p(1 -p), where p is the local frequency of allele A). They also provide some insights regarding the cause of the surprising loss of polymophism at the major-effect locus with strong selection. The simulation procedure is detailed as follows.

Populations and habitats. The species is split in two subpopulations living in two different habitats, with local carrying capacity K = 10 4 . In each habitat, individuals experience selection toward a local trait optimum θ i = (-1) i (for habitat i). Initially, the two subpopulations are at 4 5 of the local carrying capacity. The genetic information of the individuals of the initial population is set as follows. In each subpopulation, all the individuals have, at the major-effect locus, the allele whose effect is the closest to the optimum of the habitat they are in (η in habitat 2 and -η in habitat 1). The polygenic background is then set randomly and uniformly.

Genetic architecture.

We consider L = 200 unlinked loci constituting the polygenic background.

At each of these loci, two alleles segregate, having an additive effect on the trait of the individual of value σ LE √ L orσ LE √ L , where σ 2 LE is the variance at linkage equilibrum of the quantitative background. No mutation occurs at those loci. We set the variance at linkage equilibrium to σ LE = 0.1 small, so that our analysis in a small variance regime is a good approximation. (In Appendix G, we consider the same framework with a smaller number of loci involved in the quantitative background L = 50).

There is an additional locus of interest, which carries the major-effect alleles +η or -η. This locus is also unlinked to all the others and no mutation occurs at this site. Note also that the trait range, given by [-η -

σ LE √ L, η + σ LE √ L] = [-η - √ 2, η + √ 2]
extends beyond the local optima (-1,1), even in the absence of major-effects.

Life cycle.

The life cycle involves overlapping generations of small time length ∆t = 0.1. The life cycle proceeds as follows:

1. reproduction: each individual of the metapopulation chooses at random one mate within its subpopulation, and their mating produces an offspring with probability ∆t.

2. selection-competition: each individual (including offspring generated in the previous step) faces a selection-competition trial according to its trait ζ and habitat i in which they are currently living. They survive with probability exp -g∆t(ζ -θ i ) 2 -∆t Ni K and are removed otherwise (here N i denotes the size of the subpopulation i after reproduction).

3. migration: at each migration event, within each subpopulation i, a number of migrants is drawn, according to a Poisson law with parameter m∆t N i (with a hard cap of N i , which is the number of individuals currently in the subpopulation after the selection-competition step). Migrants are uniformly sampled accordingly within the subpopulation and are moved to the other deme. We stress that a given value of the migration rate m = 0.8 means that, on average, a fraction m∆t = 0.08 of the population will change deme at each generation.

Each simulation repeats this life cycle, first without migration for 100 generations of burn in (10 time units) and next with migration for N gen = 10 4 generations (10 3 time units). We model two types of initial events when migration starts: either nothing happens, and the initial state is symmetrical, or we model a sudden catastrophic loss of population in only one of the habitat during the first generation with migration, so that the initial state is asymmetrical (results shown in Appendix G). Precisely, we change the mortality of the uniform competition term only in the habitat 1, by replacing exp -∆t N1 K by exp -N1 K (which is consistent with the interpretation that this catastrophic loss of population is very abrupt). This leads to asymmetrical initial subpopulation sizes.

Qualitative results of the IBS on the stability of polymorphism with regard to increasing selection. The solid lines in the sufigures of Fig. 5 (symmetrical initial population sizes) and Fig. 9 (asymmetrical initial populations sizes) represent the median trajectories of the variance at the major-effect locus (p(1 -p), where p is the local frequency of allele A) in each habitat (gold lines for habitat 1 and navy ones for habitat 2). When the variance p(1 -p) is positive, the A/a polymorphism is maintained. In both Fig. 5 and Fig. 9, selection increases from top to bottom and the polygenic background is present in the left panel and absent in the right one. When there exists a polygenic component contributing to the trait, polymorphism at the major-effect locus is lost after some time with weak selection (g = 0.1, Fig. 5a and Fig. 9a), is maintained with intermediate selection (g = 0.5, Fig. 5c and Fig. 9c) and lost again even more quickly with strong selection (g = 1, Fig. 5e and Fig. 9e). This is qualitatively consistent with the analytical predictions displayed in Fig. 3a, where the red crosses indicate the three selection-migration set of parameters chosen for the IBS. Moreover, this phenomenon is robust with regard to initial conditions (Fig. 5 and Fig. 9), although the loss of polymophism at the major-effect locus at weak and strong selection is faster when subpopulations sizes are initially asymmetrical (Fig. 9).

Control case without polygenic background.

To confirm that the loss of polymorphism at weak and strong selections is due to the polygenic background and not to genetic drift (although drift is unlikely to have an effect under this time range of 10 3 time units with a population of order 10 4 ), we additionally run an equal number of replicates for each set of parameters without any polygenic background (L = 0, σ LE = 0). Only the major-effect alleles segregates, and this corresponds to the one-locus haploid model. Results shown in the right panel of Fig. 5 and Fig. 9 are consistent with the one-locus haploid model analysis, which states that the polymorphism at this major-effect locus is stable at all level of selection (the variance at the major-effect locus remains positive and stable).

Explanation behind the loss of polymorphism with strong and weak selection. The IBS allow us to gain some insights about the cause of the major locus polymorphism's collapse. We first consider strong selection. In particular, the dynamics of the subpopulations sizes and the local mean traits reveal that, at one point, stochasticity creates a small shift in the local mean traits. This shift is the same between bearers of A and a and in both patches, because the small segregational variance of the quantitative background binds the quantitative background values to be approximately the same for everyone (in the analysis, this is reflected by the change in variables (11) introducing δ a ε , δ A ε and δ ε ). Therefore, this shift, which is toward one of the local optima improves the adaptation of one subpopulation and is deleterious in the other one, which causes an asymmetry in subpopulation sizes, which is particularly pronounced when selection is strong (see the small figure embedded in Fig. 5e, built by selecting one of the two asymmetries for the sake of clarity). Because of this asymmetry, the migrants' flow is also asymmetrical and the larger population then undermines even more the small population by gene flow, which in turn raises the frequency of the major-effect allele favoured in the larger patch, which then further increases the disparity in population sizes among the two patches. This positive feedback loop creates a vortex related to the phenomenon of migrational meltdown identified in the quantitative genetic model of [START_REF] Ronce | WHEN SOURCES BECOME SINKS: MIGRA-TIONAL MELTDOWN IN HETEROGENEOUS HABITATS[END_REF], which eventually leads to the loss of one of the major-effect alleles.

One crucial feature of this explanation is the dynamics of the (varying) subpopulation sizes, which our eco-evo model allows us to track. To confirm this intuition, we conducted the same IBS, but with adjusting the birth rate to compensate for the deaths at every generation, effectively keeping both subpopulations at a constant size. With constant subpopulations sizes, the polymorphism at the major-effect locus is not lost with strong selection (see Fig. 8c in Appendix F). This truly highlights the role of the eco-evo framework in which subpopulations sizes are variables that are allowed to vary.

We next consider weak selection, where one might wonder if the loss of polymorphism at the majoreffect locus relies on the same mechanism. In this case, random fluctuations at the major-effect locus cause the quantitative trait to shift in the opposite direction, ensuring that the mean trait remains near 20 0 (the midpoint between the two optima). This is because migration is so high relative to selection that selection favours lineages that survive well in both patches. However, here, the small figure embedded in Fig. 5a suggests the role of varying subpopulations sizes in this phenomenon is not as important. This is confirmed by the fact that the loss of polymorphism also occurs in IBS where subpopulations sizes are kept at a constant level (Fig. 8a). This implies that the loss of polymorphism at the majoreffect locus with weak and with strong selection fundamentally differ. With strong selection, randomly generated asymmetries drive the system toward specialization in one patch and fixation of the major allele in that patch, whereas with weak selection, randomly generated asymmetries drive one major allele to fix and the quantitative trait to compensate in such a way that individuals are generalist with a mean trait near the midpoint between the patch optima.

Quantitative comparison of IBS with the continuous-time deterministic model (1).

We also ran deterministic numerical iterations of (1) to check the quantitative constistency of the stochastic IBS with the deterministic model (1). Two series were run, one for each type of initial condition (symmetrical or asymmetrical initial subpopulations sizes). The median trajectory obtained from these deterministic numerical resolutions of (1) for each set of parameters and initial conditions are displayed by the dashed lines in all the subfigures of Fig. 5 and Fig. 9. These deterministic trajectories are in excellent agreement with the ones obtained from the stochastic IBS and provide good approximations. We choose to distinguish the two types of initial conditions (asymmetrical or symmetrical) because the deterministic numerical resolutions are unequally sensitive to them. Indeed, since the environment is symmetrical, the symmetrical initial state is at an unstable edge between two symmetrical stable valleys for the interesting range of parameters and thus wanders for some time before choosing a valley to fall into. Therefore, we initialized the deterministic resolutions with the symmetrical initial state according to their respective IBS replicate states at 20% of the time before the median fixation time (140 times units for Fig. 5a and 40 time units for Fig. 5e). With the asymmetrical initial conditions (shown in Fig. 9 in Appendix G), this sensitivity is greatly reduced, and the deterministic resolutions are initialized according to their respective IBS replicate states when migration starts (0 time units). Furthermore, the numerical scheme for the resolution of the deterministic model (1) uses a splitting scheme to handle successively migration and the ecological dynamics internal to each habitat. For the latter, we use a discretization of the Duhamel' integral formula on time step of lengths ∆t for the asymmetrical initial state series and ∆t 4 for the symmetrical initial state series (with an accordingly increased number of time steps). The left panel is obtained with both a major-effect locus (η = 1/2) and a polygenic background of 200 loci, whereas only the major-effect locus is present in the right panel. For each subfigure, 20 replicates simulations were run per set of parameters, according to the setting explained in Section 3.2. In each subfigure, the solid line represents the median trajectory and the shaded area indicates the 0.2 and 0.8 quantiles. The dashed lines represent the median trajectories of the numerical resolutions of the deterministic model (1). This figure confirms that polymorphism of the major-effect locus is maintained only when selection is intermediate in strength (panel c) in presence of a polygenic background (left panel). The small figures embedded in each figure of the left panel represent the dynamics of the subpopulation sizes (N 1 and N 2 ). They highlight the qualitative difference between the loss of polymorphism at the major-effect locus with weak or strong selection, as fixation occurs without change in subpopulation size with weak selection (Fig. 5a) and as subpopulations sizes become asymmetrical with strong selection (Fig. 5e).

Discussion

Summary. In this work, we present a new eco-evo model for selection in a heterogeneous environment that combines a major-effect locus with a quantitative genetic background, without assuming that the latter is normally distributed. With this model, we aim to examine how the presence of a small quantitative background can disturb the polymorphism at the major-effect locus, which on its own would be favoured in the type of setting we consider. This model bridges a population genetic model (one-locus haploid model) with a quantitative genetic model recently studied in a heterogeneous environment [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF]). To do so, it introduces a new reproduction operator, inspired by the infinitesimal model, that encodes the inheritance of a major-effect and a quantitative background. The analysis goes deeper than previous studies, by formally justifying that the polygenic component of the trait is normally distributed around the major-effect allelic effects in a regime of small variance and hence justifying the Gaussian assumption made in [START_REF] Lande | The response to selection on major and minor mutations affecting a metrical trait[END_REF] and [START_REF] Chevin | Selective Sweep at a Quantitative Trait Locus in the Presence of Background Genetic Variation[END_REF]. To show this, we find new convex analysis arguments that leads to a separation of time scales, which allows us to study the stability of the polymorphism at the major-effect locus. We show that this polymorphism, which is maintained at intermediate selection, is subsequently lost when selection increases beyond a certain threshold, a phenomenon qualitatively confirmed by individual-based simulations. The separation of time scales' point of view offers the interpretation that the infinitesimal background slowly disrupts the rapidly established symmetrical polymorphism at the major-effect locus. Therefore, this phenomenon cannot be predicted by the one-locus haploid model (without the quantitative background). To our knowledge, this phenomenon has not yet been documented.

The importance of the eco-evo framework and the influence of small segregational variance. In the last section, we provided an explanation for our main biological result, which is the unexpected loss of polymorphism at the major-effect locus with both strong and weak selection. With strong selection, the explanation relies on two factors. First, the mean quantitative background is constrained to move similarly in both patches and for bearers of A and a because of the small segregational variance. This implies that any slight shift of the mean quantitative background necessarily increases local adaptation to one patch and decreases local adaptation to the other. Consequently, the latter creates an asymmetry in subpopulation sizes, one being better adapted than the other. This asymmetry is significant when local selection is strong. The larger subpopulation then sends relatively more migrants to the other patch, undermining the local adaptation there even more, which contributes to raise the frequency of the allele favoured in the now larger patch everywhere. In turn, the combination of increasing specialization and increasing disparity in population sizes (and therefore migrant production results) in a vortex that can be identified as a migrational meltdown (coined in Ronce and Kirkpatrick 2001).

Therefore, one can observe that this phenomenon specifically relies first on our eco-evo framework, which allows us to track the dynamics of subpopulation sizes. This highlights the importance of considering eco-evo dynamics when dealing with strong selection (which can heavily impact population sizes), as this loss of polymorphism at the major-effect locus with strong selection would not be captured by a more traditional approach in population genetics that considers the population sizes constant. Second, this phenomenon of loss of polymorphism at the major-effect locus also relies crucially on the small segregational variance of the quantitative background, which is linked to the very small effect sizes of sufficiently many alleles. It is indeed worth noting that, if the segregational variance of the quantitative background is relatively large, then the first step of our explanation ("the mean quantitative background is constrained to move similarly in both patches and between bearers of A and a") does not hold and the mean quantitative background can shift in opposite directions in the two patches, improving local adaptation in both patches. The impact of bimodality in the quantitative trait, with mean trait values in each patch near the local optimum, on the stability of polymorphism at the major-effects locus deserves further attention.

Robustness.

To assess the robustness of the swamping phenomenon that we identified, we performed various individual-based simulations. We found that these are in excellent quantitative agreement with our analysis. They also connect our framework to the evolution of an explicit genetic architecture, which provides a practical translation of the small variance regime that underlies our study. This is important, because we have only shown that our results hold in this small variance regime. In particular, they might be different under parameter ranges that violate this regime, for example under low or no migration (meaning, at a level of comparable order as the small variance). Moreover, because the trajectories of the individual-based simulations are consistent with the deterministic trajectories produced by our model, we can validate essential assumptions underlying our model, mostly the constancy of the small segregational variance for the quantitative background. The latter requires enough loci (L) with relatively small effects (± σLE √ L ), so that the segregational variance of the quantitative background (lower than σ 2 LE ) remains small while the phenotypic range produced by the polygenic background alone -σ LE √ L, σ LE √ L spans well beyond the local optimal traits. The last condition is necessary to ensure genotypic redundancy (see also [START_REF] Yeaman | Evolution of polygenic traits under global vs local adaptation[END_REF], so that well adapted mates with similar phenotypes have on average relatively different genotypes, which in turn ensures that the variance of their offspring does not depend too much on their traits. In our simulations presented in Section 3.2, we showed that L = 200 and σ LE = 0.1 produced very similar trajectories to our deterministic model. Furthermore, in Appendix G, we even lowered the number of loci to L = 50 and increased the segregational variance parameter σ LE = 0.2 to assess the robustness of our conclusions with regard to less favorable parameters, with the same conclusions.

Complete analytical outcomes.

The analysis performed in Section 3.1 is centered on the persistence of polymorphism at the major-effect locus. As stated in Remark 1, the loss of polymorphism by fixation would lead to the dynamics of the quantitative background alone, as covered in Dekens (2022). Hence, Fig. 6 complements Fig. 3 (for η = 0.5 and varying migration and selection). First, it displays the region of parameters where where polymorphism at the major-effect locus would be maintained, leading to a bimodal trait distribution in the metapopulation (in yellow, corresponding to the analogous region in Fig. 3a). For the rest of the parameters, either the population goes extinct (purple region, when both migration and selection are too strong), or one major-effect allele fixes and the other one is lost. Therefore, in this case, the equilibria are given by the analysis done in Dekens (2022) (as anticipated by Remark 1). These equilibria are monomorphic and are of two types. First, for bounded selection, there exists a critical threshold in the migration rate under which the polymorphism at the major-effect locus is lost due to the strong blending effect of migration. In that case of strong relative migration, the metapopulation occupies equally the two habitats as generalist and its trait distribution concentrates around the midpoint between the two habitats' optima. It is therefore constituted by individuals that do not suffer from a loss in fitness when migrating between the patches. This corresponds to the symmetrical monomorphic equilibrium (see [START_REF] Ronce | WHEN SOURCES BECOME SINKS: MIGRA-TIONAL MELTDOWN IN HETEROGENEOUS HABITATS[END_REF][START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF], where the population can be qualified as generalist (green region). Second, for bounded migration rates, the polymorphism at the major-effect locus lost above a certain threshold in selection. In this case of relative strong selection (blue region), the metapopulation becomes specialized in one of the two patches, and deserts the other (at the exception of a few migrants who are too maladapted to establish there). This type asymmetrical equilibrium, highlighted as a source-sink scenario in [START_REF] Ronce | WHEN SOURCES BECOME SINKS: MIGRA-TIONAL MELTDOWN IN HETEROGENEOUS HABITATS[END_REF], was analytically derived in [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF].

Filling a methodological gap.

In population genetics, one-locus or two-locus models in heterogeneous environments have been well studied (Nagylaki and Lou 2001; Bürger and Akerman 2011), with a nuanced picture when including the effect of drift (Yeaman and Otto 2011). A two-deme two-locus model is analysed in Geroldinger and Bürger (2014), which in particular shows that a concentrated genetic architecture (a major-effect locus and a tightly linked minor one) maintains polymorphism (full or single-locus) even under high migration rates when selection acts in opposite directions in the two patches. Increasing the number of loci quickly leads to analytical complexity too great for general study. There also exist multi-loci models in heterogeneous environments [START_REF] Lythgoe | Consequences of gene flow in spatially structured populations[END_REF][START_REF] Szép | Polygenic local adaptation in metapopulations: A stochastic eco-evolutionary model[END_REF]), but they focus on equal allelic effects. On the other end of the spectrum, quantitative genetic models do not typically account for additional discrete major-effect allelic effects on the focal quantitative trait (for sexually reproducing populations in heterogeneous environment, see Ronce and 2022) (equivalent to this model upon loss of polymorphism). For bounded selection, when migration increases, there is a threshold over which the polymorphism at the major-effect locus is lost due to the blending effect of migration (consistent with Yeaman and Whitlock 2011). The population then becomes equally maladapted to both habitats (generalist -symmetrical monomorphism, in the green region). For this specific major-effect allelic effect η = 0.5, there exists additionally an interval of selection strength (≈ [0.7, 1]) for which the major polymorphism might not be stable at all migration rates below the critical threshold. This phenomenon does not seem to hold when the major effect is larger (see Fig. 3b). For bounded migration (below the threshold rate over which the strong migration blending hampers the major polymorphism), when selection strength increases, the polymorphism at the major-effect locus (yellow region) is lost, and the population becomes adapted to one of the two habitats (specialist -asymmetrical monomorphism, in the blue region). As this figure is obtained in the small segregational variance regime (which should remain smaller than the other parameters of the system for the analysis to be valid), we warn that the outcomes displayed in the vicinity of the x-axis (very small migration rates) might not correspond to the limit when the migration rate is 0 (no migration).

To our knowledge, the first model that bridges this gap between quantitative traits and discrete loci appears in [START_REF] Lande | The response to selection on major and minor mutations affecting a metrical trait[END_REF]. In this work, the author considers the dynamics of a single major-effect locus where two alleles segregate along with a polygenic background, in a diploid panmictic population subjected to a sudden change of environment. He models the influence of the polygenic background on the trait by assuming that, among bearers of the same major-effect allele, the trait distribution is Gaussian, centered around the effect of the the major-effect allele on the trait. This study opened the way for more recent work on the genetic architecture of adaptation of panmictic populations in a suddenly changing environment, where the central question is whether this adaptation is due to majoreffect allelic sweeps or to subtle shifts in the frequency of many small effect alleles. In [START_REF] Chevin | Selective Sweep at a Quantitative Trait Locus in the Presence of Background Genetic Variation[END_REF], the authors extend the framework of [START_REF] Lande | The response to selection on major and minor mutations affecting a metrical trait[END_REF] to include less specific selection functions than exponential ones. Subsequent studies (Vladar and Barton 2014; Jain and Stephan 2017) explicitly model the short-term dynamics of a polygenic trait at mutation-selection balance, following a sudden change of environment. They show that there exists a sharp threshold in allelic sizes below which polymorphism remains and above which fixation occurs. Lately, in a similar context, Höllinger, Pennings, and Hermisson (2019) propose an extension to take genetic drift into account on the dynamics of adaptation with a polygenic binary trait under mutation-selection balance. However, all those works from Lande (1983) to Höllinger, Pennings, and Hermisson (2019) study panmicitic populations, without spatial structure, even though spatial heterogeneities are known to generate gene flow, which indirectly shapes genetic architecture through local adaptation (see [START_REF] Yeaman | THE GENETIC ARCHITECTURE OF ADAP-TATION UNDER MIGRATION-SELECTION BALANCE: THE GENETIC AR-CHITECTURE OF LOCAL ADAPTATION[END_REF], or below for more details). Moreover, they focus solely on the dynamics of the allelic frequencies without considering their coupling with population size dynamics, assuming it to be constant.

In this paper, we presented a composite framework between population and quantitative genetics aiming at going beyond these methodological limitations. Our model and methodology allows us to study analytically the eco-evo dynamics of a sexually reproducing population characterized by a composite trait resulting from the interplay between a few major-effect loci and a quantitative polygenic background, in spatially heterogeneous environments (migration-selection balance). We want to emphasize that, by "eco-evo dynamics", we mean that we study both the ecological and evolutionary dynamics of the local trait distributions and therefore do not assume that the sizes of the populations remain constant; rather, they are variables of the system. This modelling choice is crucial, because the migrational meltdown phenomenon provoking the loss of polymorphism at the major-effect locus with strong selection relies on the building of asymmetrical subpopulation sizes.

The role of the Gaussian assumption of quantitative trait values.

In our work, we justify the Gaussian assumption made by [START_REF] Lande | The response to selection on major and minor mutations affecting a metrical trait[END_REF] and [START_REF] Chevin | Selective Sweep at a Quantitative Trait Locus in the Presence of Background Genetic Variation[END_REF] to model the background polygenic effect on the trait via a framework that does not make a priori assumptions on the within-population distribution. Instead, our model relies on an extension of the standard infinitesimal model [START_REF] Fisher | The Correlation between Relatives on the Supposition of Mendelian Inheritance[END_REF]) that encodes both the inheritance of the quantitative background and the major-effect alleles. Analytical progress is possible in a regime of small segregational variance for the quantitative component of the trait, despite not specifying the shape of the trait distribution. It relies on the fact that the variance introduced at each event of reproduction by the quantitative background is small compared to the discrete allelic effects at the major-effect locus. This allows us to use a methodology developed by [START_REF] Diekmann | The dynamics of adaptation: An illuminating example and a Hamilton-Jacobi approach[END_REF], meant to study trait distributions concentrated as Dirac masses, to justify that assuming Gaussian distributions of quantitative trait values is valid (Section 2.1). Moreover, this Gaussian approximation appears here sufficient to capture the phenomenon of migration meltdown with strong selection that we identify through the rest of our analysis, as higher order moments do not seem to have a significant influence on it.

Extensions to more complex population genetic models. The model and the line of methods that we use in this paper are quite robust. We thus provide a comprehensive toolbox at the end of this manuscript (Appendix A), to describe how to apply the method more broadly. In particular, the toolbox is meant to indicate how to extend the method to more complex population genetic models by adding a quantitative background. It relies on Proposition B.1, which justifies that carrying the analysis under the Gaussian assumption of quantitative trait distribution is valid (in the regime of small variance indicated). In Appendix A, we detail the hypotheses that the population genetic models must satisfy in order to use it (see Appendix B for details and examples).

Further prospects. The loss of the polymorphism at the major-effect locus with strong selection in a symmetrical heterogeneous environment, where one might think that it is most favoured, illustrates the value of our method. However, two natural questions stem from our work:

1. Would the stability region of polymorphism at the major-effect locus shrink as much in the presence of a quantitative background when considering asymmetrical levels of selection/migration between the two patches? Our analysis suggest that it should, and this can be investigated through an extension of the last step of our analysis.

2. Would this phenomenon hold if mutations can accumulate at the major-effect locus? Figure 3 for example suggests that polymorphism at the major-effect locus would persist over a wider range of parameters if the alleles at the major-effect locus evolve to match the difference in optima. This possibility was indicated by the numerical findings of Yeaman and Whitlock (2011), who found the emergence of tightly linked clusters of major-effect loci underlying local adaptation for intermediate migration rates.

A Toolbox: How to study the interplay between a quantitative background and a finite number of major-effect loci.

The aim is to study the interplay between a quantitative background and a finite number of major-effect loci. We start with a population genetic model. Let us consider K different genotypes A (k) which have genotypic effects on the phenotype a (k) (we use the index k to indicate genotypes). For our method to be applied, the genotypes should verify two hypotheses H1 and H2 described in Appendix B. The metapopulation lives in a heterogeneous environement of I patches (we use the index i to indicate location). We denote the population of patch i carrying genotype k by N (k) i . Let us denote the system of equations that describes the dynamics of the genotypic local population sizes: d N dT = Gā N (T ) and of a viable stable equilibrium N * . We recall that N * is an equilibrium of the system if Gā ( N * ). This equilibrium is viable if all the population sizes are non-negative, and at least one is positive. Its local stability is determined by standard linear analysis (sign of the real parts of the eigenvalues of the system's Jacobian).

Let us modify the previous population genetic framework to include the effect of a quantitative background on the trait, generically denoted z. While previously, all individuals carrying the same genotype A (k) shared the same phenotype, now their phenotypes can differ due to the quantitative background they present. Consequently, among individuals of the same patch k carrying the same major genotype A (k) , we distinguish those sharing the same quantitative background z, and denote their number n (k) i (z):

A (k) ⇝ (A (k) , z) a (k) ⇝ a (k) + z N (k) i ⇝ n (k) i (z).
The PDE system that we obtain on the trait distributions n (k) i is not easily analysed. That is why we provide a five steps plan in order to guide the analysis when the diversity introduced by the segregation of the quantitative component of the trait is small compared to the variance generated by the majoreffect loci (H6 -regime of small variance):

1. First, we operate a scaling of time according to the regime of small variance. It anticipates on the separation of time scales such that the major-effect allelic frequencies change rapidly, followed by the slow changes of the quantitative components (see step 3).

2. In this regime of small variance, we can justify the Gaussian approximation of the local genotypic distributions n (k) i centered at the same mean and the same variance ε 2 , thanks to Proposition B.1, as soon as the assumptions (H1) and (H2) are satisfied (see Appendix B) and every genotypic population randomly mates with themselves and every other genotypic population (H3) (the latter excludes for example models that differentiate sexes). This guides the intuition toward which change of variables to perform in order to get a system separating time scales explicitly (see Step 3). We emphasize that Proposition B.1 is crucial to be able to apply this method.

3. From the PDE system on the distributions n (k) i , we can deduce the ODE system of their moments. Since we have justified the Gaussian approximation for all local genotypic distributions n (k) i , the new system is closed in the regime of small variance ε 2 ≪ 1, and only involves the dynamics of the genotypic local sizes of populations N (k) i and the genotypic local mean quantitative components z (k) i . 4. This step aims at obtaining a system that explicitly separates time scales, in order to ultimately reduce the complexity of the analysis. It a technical change of variables, which is guided by the formal analysis of the step 1 (mean quantitative components roughly the same within patches), and the intuition that migration has a strong blending effect between patches in the small variance regime (which would result in the mean quantitative components roughly being equal between patches). These considerations bring the following new variables replacing the genotypic local mean quantitative component z

(k) i : ⋄ for each genotype 1 ≤ k ≤ K, δ (k)
i,ε is the difference in the mean quantitative component of the genotypic population k between the patch i + 1 and patch i (1 ≤ i ≤ I -1). Dividing by ε 2 comes from the intuition given above.

⋄ for each genotype 1 ≤ k ≤ K -1, δ (k) ε
is the difference between the mean quantitative component averaged across patches of genotypic population k + 1 and k. Dividing by ε 2 comes from the intuition given in Step 1.

⋄ Z ε is the overall mean quantitative component across patches and major genotypes. It is the slow evolving variable.

Rewriting the dynamics of the genotypic local population sizes Nε along these new variables δε and Z ε delivers a system in which all the differential equations are multiplied by ε 2 (fast dynamics of Nε and δε ) except the one governing the dynamics of Z ε (slow dynamics).

To finally complete the separation of time scales and obtain the limit system by letting ε 2 vanish, it is sufficient to show that at each value Z of the slow variable, the fast time-scale equilibria N , δ are stable, for example by using the Routh Hurwitz criterion for linear analysis on the Jacobian Jac Gā N , δ .

5. The last step to determine the stability of the global equilibria of the full system of the genotypic population with the influence of the quantitative background N * , δ * , Z * , consists in applying the formula given in the last box (see the next page).

Toolbox: How to study the interplay between a quantitative background and a finite number of major-effect loci dynamics.

The stadium:

I patches P i (1 ≤ i ≤ I)
The teams:

K different genotypes A (k) (1 ≤ k ≤ K)
Vector of genotypic effects on phenotype: ā = (a (1) , ..., a (k) ) Matrix of local genotypic population sizes: N = N (k) i

Pop. gen. model:

d N dT Gā N (T )
Pop gen. analysis:

(i) Viable equilibria: G ā N * = 0 and N * > 0.

(ii) Stability: eigenvalues of Jac Gā N * in open left plane.

Population genetic model

The new players:

(i) Quantitative background z (ii) Individuals carrying genotype i and a quantitative background z have a phenotype z + a (k) .

(iii) Distribution in patch k : n

Work hypotheses:

H1 -H2 (reflexivity and irreducible graph -see App.B)

H3 every genotypic population reproduces randomly with themselves and every other in the same patch H4 inheritance of the quantitative background in accordance with the infinitesimal model with segregational variance σ 2 .

H5 the quantitative background is unlinked to

A (k)
H6 σ 2 ≪ min a (k) 2 : small variance regime.

Composite model combining population and quantitative genetics

(i) n (k) i,ε (z) ≈ N (k) i,ε × Gauss z (k) i,ε , ε 2 (ii) z (k) i,ε ≈ z (l) i,ε 2) ODE system of moments (z ε := (z (k) i,ε )):    ε 2 d Nε dt = G ā Nε (t), zε (t) , ε 2 dzε dt = F ā Nε (t), zε (t) .
3) Slow-fast analysis:

(i) Change in variables: δ (k) i,ε = z (k) i+1,ε -z (k) i,ε 2ε 2 [K(I -1)] ; δ (k) ε = i z (k+1) i,ε -z (k) i,ε 2Iε 2 [(K -1)] ; Z ε = k,i z (k) i,ε

K×I

(ii) Slow-fast system:

   ε 2 d[ Nε, δε] dt = G ā Nε (t), δε (t), Z ε , dZε dt = F ā Z ε , Nε , δε .
(iii) Separation of time scales (via stability of zeros of G ā by Routh-Hurwitz criterion on Jac Gā ( N , δ))

   0 = G ā N , δ, Z , dZ dt = F ā Z, N , δ .
4) Analysis of the limit system:

(i) Viable equilibria: G ā N * , δ * , Z * = F ā Z * , N * , δ * = 0, N * > 0 (ii) Stability: ∇ N , δ F ā • Jac Gā N , δ -1 ∂ Z G ā Z * , N * , δ * > 0.
Steps to apply the analysis B Generalization of Proposition 2.1 for more complex genotypes.

To state a generalization of Proposition 2.1, we first need to specify the targeted scope of population genetic models. Let us consider K different genotypes A (k) that satisfies the following hypotheses relating to how they interact with each other regarding the genotypes of their offspring:

H1 Reflexivity: For all k ∈ (1, K), the offspring of two parents with same genotype A (k) has a positive probability to be of genotype A (k) .

H1 is a natural hypothesis when considering either haploid or diploid populations, even with non-Mendelian processes (genetic linkage/recombination, gene drives), provided that they are not too extreme (lowering the probability of inheriting a certain genotype is fine as long as it does not cancel it).

The second hypothesis is more conveniently apprehendable by considering the graph G whose nodes are the genotypes A (k) . A vertex links two nodes A (k) and A (l) if and only if there exists a positive probability that their offspring has genotype A (k) or A (l) .

H2 Irreducible graph: For all (k, l) ∈ (1, K) 2 , there exists a path of vertices of G connecting A (k) and A (l) .

This last hypothesis is satisfied by any haploid models, regardless of how many loci are considered, because an offspring can inherit all their alleles from only one parent. Consequently, in that case, every node of the graph is connected to every other. In diploid models, where an offspring can have a different genotype from both its parents, which vertices of the graph G exist is not clear. However, for example, we can show that the graph corresponding to a diploid model, with L loci and two alleles at each loci, is connected according to H2. Indeed, each genotype is directly connected to any other that differs from it from just one allele at one locus. Nevertheless, the interest of H2 is that it is very easy to verify whether it is satisfied given any particular model.

To state the proposition that generalizes (Proposition 2.1), we first need to define the index set of couples that can yield an offspring with a particular genotype. For k ≤ K, we denote it by C (k) , where (l, k) ∈ C (k) if and only if parents with genotypes A (l) and A (k) can produce an offspring with genotype A (k) . The following proposition characterizes the genotypic functions u A (k) that respect the following constraints analogous to C ∀k ≤ K, ∀z ∈ R, max

(l,k)∈C (k) sup z1,z2 u A (k) (z) -z - z 1 + z 2 2 2 -u A (l) (z 1 ) -u A (k) (z 2 ) = 0. (C')
Proposition B.1. Suppose that H1 and H2 are satisfied. For k ≤ K, we consider u A (k) a real valued non-negative function whose zero set is non-empty and of measure 0 (for example, is finite). If {u A (k) , k ≤ K} respects (C'), then there exists z * ∈ R such that for all k ≤ K:

∀z ∈ R, u A (k) (z) = (z -z * ) 2 2 . Proof. 0) u A (k)
is continuous and has right and left derivatives everywhere. For k ≤ K and z ∈ R, we have:

u A (k) (z) -z 2 = min (l,k)∈C (k) inf z1,z2 -z(z 1 + z 2 ) + z 1 + z 2 2 2 + u A (l) (z 1 ) + u A (k) (z 2 ) . ( 16
)
Therefore, u A (k) (z) -z 2 is concave as infimum of affine functions, and thus continuous and has right and left derivatives. 1) u A (k) cancels only once. Let us fix k ≤ K. Let us suppose that u A (k) has two zeros z * 1 ̸ = z * 2 . H1 implies that (k, k) ∈ C (k) . Then, we deduce from (C') that:

u A (k) (z) ≤ inf z1,z2 z - z 1 + z 2 2 2 + A (k) (z 1 ) + u A (k) (z 2 ).
In particular, for z =

z * 1 +z * 2 2 , z 1 = z * 1 , z 2 = z * 2 , we obtain u A (k) ( z * 1 + z * 2 2 ) ≤ 0. ( 17
)
As u A (k) is non-negative, the midpoint between two zeros of u A (k) is also a zero of u A (k) . u A (k) is also continuous from the previous point, therefore, we deduce that u

A (k) cancels on [z * 1 , z * 2 ]
. The latter violates the assumption that u A (k) has a zero set of measure 0. Because it is also not empty, we get that u A (k) cancels exactly once, in a point that we denote z * k .

2) The zero of u A (k) coincides with the zero of u A (l) : z * k = z * l . First, let us consider the case where (k, l) ∈ (1, K) 2 is such that A (k) and u A (l) are linked by a vertex in the graph G. Then, we deduce that (k, l) ∈ C (k) or (k, l) ∈ C (l) . We can assume the first without loss of generality. Similarly as the first part of the proof, we deduce that

u A (k) (z) ≤ inf z1,z2 z - z 1 + z 2 2 2 + u A (k) (z 1 ) + u A (l) (z 2 ).
Consequently, the midpoint between z * k and z * l is a zero of u A (k) , which is necessarily z * k , which implies that z * k = z * l . Let us now show the same for every couple (k, l) not necessarily linked by a vertex in G. H2 implies that there exists a path of vertices between u A (k) and u A (l) . As we showed that for every pair of nodes connected by a vertex, the zeros of their function is the same point, that property also holds for the extremities of the path of vertices, hence z * k = z * l . We denote z * the common zero.

3) Convex Legendre conjugates ûA (k) (y) = sup z (z -z * )y -u A (k) (z). Let us show that (C')

implies that the convex Legendre conjugate satisfies

∀y ∈ R, ûA (k) (y) = y 2 4 + max (l,k)∈C (k) ûA (l) y 2 + ûA (k) y 2 . ( 18
)
Using ( 16) and commuting the sup, we obtain, for y ∈ R,

ûA (k) (y) = sup z (z -z * )y -min (l,k)∈C (k) inf z1,z2 z - z 1 + z 2 2 2 + u A (l) (z 1 ) + u A (k) (z 2 ) = max (l,k)∈C (k) sup z1,z2 -u A (l) (z 1 ) -u A (k) (z 2 ) + sup z (z -z * )y -z - z 1 + z 2 2 2 . ( 19
)
Moreover, a straight-forward calculus shows that the sup

z is reached at z = y+z1+z2 2 , which leads to sup z (z -z * )y -z - z 1 + z 2 2 2 = y + z 1 + z 2 2 -z * y - y 2 4 = y 2 4 + (z 1 -z * ) y 2 + (z 2 -z * ) y 2 . ( 20
)
Combining ( 20) and (19) (the fact that z A = z a = z * plays a crucial part for the crossed term) leads to (18). Moreover, we obtain classically that:

ûA (k) (y) ≥ (z * -z * )y -u A (k) (z * ) = 0 = ûA (k) (0) (21)
.

4) max

k≤K ûA (k) : y → y 2 2 . We obtain (18) that:

∀y ∈ R, max k≤K ûA (k) (y) = y 2 4 + max k≤K max (l,k)∈C (k) ûA (l) y 2 + ûA (k) y 2 . ( 22
)
For y ∈ R, let k 0 ≤ K be such that max k≤K ûA (k) y 2 = ûA k 0 y 2 . H1 implies in particular (k 0 , k 0 ) ∈ C (k0) and therefore, the maximum of the right-hand side of ( 22) is reached in 2 ûA k 0 y 2 . Consequently, we deduce that

∀y ∈ R, max k≤K ûA (k) (y) = y 2 4 + 2max k≤K ûA (k) y 2 . ( 23
)
Moreover, one can notice that

∀y ∈ R, max k≤K ûA (k) (y) = max k≤K max z∈R (z -z * )y -u A (k) (z) = max z∈R (z -z * )y -min k≤K u A (k) (z) = ˆ min k≤K u A (k) (y).
Therefore, max k≤K ûA (k) is a convex continuous function that has left and right derivative everywhere, in particular in 0. Hence, iterating (23) implies first that:

∀y > 0 (resp. < 0), max k≤K ûA (k) (y) = y 2 2 + β y (resp. α y), ( 24
)
where (α, β) = max k≤K ûA (k) ′ (0 -), max k≤K ûA (k) ′ (0 + ) . From (21), we deduce that the α ≤ 0 ≤ β. Since max k≤K ûA (k) is the convex conjugate of min k≤K u A (k) , we compute that the convex bi-conjugate of min

k≤K u A (k) is z →      (z-z * -α) 2 2 if z < z * + α 0 if z * + α ≤ z ≤ z * + β (z-z * -β) 2 2 if z > z * + β. ( 25
)
As the convex bi-conjugate of min k≤K u A (k) is the lower convex envelope of min k≤K u A (k) , the two of them are equal at the extremal points of its graph, namely for z = z * + α and z = z * + β. We deduce from (25) that min

k≤K u A (k) (z * + α) = min k≤K u A (k) (z * + β) = 0.
Since all the u A (k) , k ≤ K only cancels for z = z * , we obtain that α = β = 0 and (24) yields that max k≤K ûA (k) : y → y 2 2 .

5) max

k≤K ûA (k) = min k≤K ûA (k) . First let us state that min k≤K ûA (k) is continuous as minimum of a finite number of continuous functions and that it is non-negative and reaches its minimum in 0, with min k≤K ûA (k) (0) = 0 (from (21)). Moreover, (18) implies that

∀y ∈ R, min k≤K ûA (k) (y) y 2 4 + 2 max k≤K ûA (k) y 2 = y 2 2 . ( 26
)
Therefore min k≤K ûA (k) has left and right derivatives in 0, and min

k≤K ûA (k) ′ (0 + ) = min k≤K ûA (k) ′ (0 -) = 0.
Furthermore, (18) also implies that

∀y ∈ R, min k≤K ûA (k) (y) ≥ y 2 4 + 2 min k≤K ûA (k) y 2 .
Iterating the last inequality, and knowing that min

k≤K ûA (k) (0) = min k≤K ûA (k) ′ 0 + = min k≤K ûA (k) ′ 0 -= 0, leads to ∀y ∈ R, min k≤K ûA (k) (y) ≥ y 2 2 = max k≤K ûA (k) (y).
Consequently, we deduce that min k≤K ûA (k) = max k≤K ûA (k) .

End of proof.

The last result implies that

∀k ≤ K, ∀y ∈ R, ûA (k) (y) = max k≤K ûA (k) (y) = y 2 2 .
From the latter we compute the bi-conjugates ûA (k) : z → (z-z * ) 2

2

. Since z → z-z * 2 is strictly convex and it is the lower convex envelope of u A (k) , we obtain that

∀k ≤ K, ∀z ∈ R, u A (k) (z) = (z -z * ) 2 2 .
C Formal justification of the constraints (C) on the main terms u A 0 and u a 0 We drop the index i indicating the habitat and the time dependence t for this appendix for the sake of simpler notations. Let us first formally justify that U A 0 and U a 0 are positive almost everywhere and cancelling somewhere. As we are interested in the maintenance of the polymorphism at the major-effect locus, we consider that no major-effect allele has yet fixed. Hence, N A ε and N a ε need to remain positive and bounded when ε vanishes. Using the Hopf-Cole transforms on n A ε and n a ε (5) along with the formal Taylor expansions (6) on U A ε and U a ε leads to

N A ε = R n A ε (z ′ ) dz ′ = R 1 √ 2πε e -U A ε (z ′ ) ε 2 dz ′ = R 1 √ 2πε e - u A 0 (z ′ ) ε 2 e -u A 1 +ε 2 v A ε dz ′ . ( 27
)
If we assume that the residues u A 1 and v A ε stay bounded when ε vanishes (as Calvez, Garnier, and Patout (2019) suggests it), then (27) implies that u A 0 must be non-negative for N A ε to remain bounded when ε vanishes. For N A ε not to vanish asymptotically, u A 0 must cancel. Moreover, for any interval I ⊂ R, u A 0 cannot cancel on I, or we would have:

N A ε > I 1 √ 2πε e 1 ε 2 e -u A 1 +ε 2 v A ε dz ′ → +∞.
So u A 0 is positive almost everywhere, and cancelling The same holds for u a 0 . Now, for determining the constraints (C), let us notice that if we divide the right-hand side of the first equality of (3) by n A ε (z), the reproduction term

B A ε (n A ε ,n a ε )(z) n A ε (z)
has to remain positive and bounded for all z ∈ R when ε vanishes for the effect of reproduction to remain well-balanced with selection, migration and competition. We assume henceforth that (6) is the correct ansatz (as suggested by Calvez, Garnier, and Patout 2019). Using the Hopf-Cole transforms on n A ε and n a ε (5) along with the formal Taylor expansions (6) on U A ε and U a ε in (4) leads to

B A ε (n A ε , n a ε )(t, z) n A ε (z) = B A ε (n A ε , n a ε )(z) 1 √ 2πε e - u A 0 (z) ε 2 e -u A 1 (z)+O(ε 2 ) = √ 2 N A ε × R 2 exp 1 ε 2 u A 0 (z) -z - z1 + z2 2 2 -u A 0 (z1) -u A 0 (z2) exp u A 1 (z) -u A 1 (z1) -u A 1 (z2) + O(ε 2 ) dz1dz2 + R 2 exp 1 ε 2 u A 0 (z) -z - z1 + z2 2 2 -u A 0 (z1) -u a 0 (z2) exp u A 1 (z) -u A 1 (z1) -u a 1 (z2) + O(ε 2 ) dz1dz2 .
As N A ε remains bounded and does not vanish asymptotically, we need the maximum of the two integrals nor to vanish, nor to diverge to infinity when ε vanishes, for all z ∈ R. Therefore the maximum of the terms into brackets that are multiplied by 1 ε 2 needs to be null for all z ∈ R:

∀z ∈ R, max sup z1,z2 u A 0 (z) -z - z 1 + z 2 2 2 -u A 0 (z 1 ) -u A 0 (z 2 ), sup z1,z2 u A 0 (z) -z - z 1 + z 2 2 2 -u A 0 (z 1 ) -u a 0 (z 2 ) = 0,
which is the first constraint of (C). The same holds for

B ⊣ ε (n a ε ,n A ε )(z) n a ε (z)
, which gives the second constraint of (C).

D Slow-fast analysis underlying the separation of time scales

Under the change of variable (11), the system (10) is equivalent to the following:

                                                                                            ε 2 d N a ε,i dt = λ i-1 N a ε,i -N A ε,i + N a ε,i N a ε,i -g i Z ε + η a -(-1) i 2 N a ε,i +α (-1) j m j N a ε,j -m i N a ε,i + O(ε 2 ), ε 2 d N A ε,i dt = λ i-1 N A ε,i -N A ε,i + N a ε,i N A ε,i -g i Z ε + η A -(-1) i 2 N A ε,i +α (-1) j m j N A ε,j -m i N A ε,i + O(ε 2 ), ε 2 d δ a ε dt = g 1 + g 2 + (g 1 -g 2 ) (Z ε + η a ) + δε 2 N A ε,2 N a ε,2 +N A ε,2 - N A ε,1 N a ε,1 +N A ε,1 -δ a ε m2α N a ε,2 N a ε,1 + m1N a ε,1 α N a ε,2 + δ A ε -δ a ε 4 N A ε,2 N a ε,2 +N A ε,1 + N A ε,1 N a ε,1 +N A ε,1 + O(ε 2 ), ε 2 d δ A ε dt = g 1 + g 2 + (g 1 -g 2 ) (Z ε + η A ) + δε 2 N a ε,1 N a ε,1 +N A ε,1 - N a ε,2 N a ε,2 +N A ε,2 -δ A ε m2αN A ε,2 N A ε,1 + m1N A ε,1 αN A ε,2 + δ a ε -δ A ε 4 N a ε,2 N a ε,2 +N A ε,2 + N a ε,1 N a ε,1 +N A ε,1 + O(ε 2 ), ε 2 d δε dt = -δε 2 -(g 1 + g 2 ) η A -η a 2 + δ A ε 2 αm2N A ε,2 N A ε,1 - m1N A ε,1 αN A ε,2 - δ a ε 2 αm2N a ε,2 N a ε,1 - m1N a ε,1 αN a ε,2 +O(ε 2 ), dZε dt = (g 2 -g 1 ) -(g 1 + g 2 ) Z ε + η A +η a 2 + δε 2 N A 1,ε N A 1,ε +N a 1,ε + N A 2,ε N A 2,ε +N a 2,ε -1 + δ a ε 2 αm2N a ε,2 N a ε,1 - m1N a ε,1 αN a ε,2 + δ A ε 2 αm2N A ε,2 N A ε,1 - m1N A ε,1 αN A ε,2 + δ A ε -δ a ε 4 N A 2,ε N A 2,ε +N a 2,ε - N A 1,ε N A 1,ε +N a 1,ε + O(ε 2 ). ( 28
)
The system (28) can be recasted more compactly into P ε . The main slow-fast analysis result is Theorem 2.1, which states the convergence of P ε towards a limit system P 0 which separates ecological and evolutionary time scales. The arguments of the proof of Theorem 2.1 are similar to the analogous theorems proved in [START_REF] Levin | Singular Perturbations of Non-linear Systems of Differential Equations and an Associated Boundary Layer Equation[END_REF] and [START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF]. The proof requires some preliminaries results, particularly of stability, to which we dedicate the rest of this section. The structure of this section is represented in Fig. 7. In the rest of this section, we first solve the slow manifold algebraic system G( Ȳ , Z) = 0, showing that there can only exist one instantaneous ecological equilibrium at a given Z ∈] -1, 1[ (Proposition D.1 and Proposition D.2). Surprisingly, this resolution is easier than the analogous one in Dekens (2022) (see Remark 5). Next, in Appendix D.1.2, we show a stability criterion of the slow manifold (Proposition D.3 and Proposition D.4), which justifies the separation of time scales approach.

D.1 Analysis of the fast equilibria.

The fast equilibria, for Z ∈] -1, 1[, are defined as the solutions Ȳ to the algebraic system G( Ȳ , Z) = 0, or equivalently seven equations that we group in two subsystems S 0 (Z) and S L (Z):

         α m 2 N a 2 -m 1 N a 1 + N a 1 1 -(N a 1 + N A 1 ) -g 1 (Z + η a + 1) 2 = 0, α m 2 N A 2 -m 1 N A 1 + N A 1 1 -(N a 1 + N A 1 ) -g 1 (Z + η A + 1) 2 = 0, m1 α N a 1 -m 2 N a 2 + N a 2 λ -(N a 2 + N A 2 ) -g 2 (Z + η a -1) 2 = 0, m1 α N A 1 -m 2 N A 2 + N A 2 λ -(N a 2 + N A 2 ) -g 2 (Z + η A -1) 2 = 0. (S 0 (Z)) J S L   δ δ A δ a   =   (g 1 + g 2 ) η A -η a 2 -(g 1 + g 2 ) + (g 2 -g 1 )(Z + η A ) -(g 1 + g 2 ) + (g 2 -g 1 )(Z + η a )   , (S L (Z)) (P ε )             
Fast dynamics (allelic subpopulations sizes and spatial discrepancies between mean infinitesimal parts), Slow dynamics (mean infinitesimal part). where:

J S L :=        -1 2 1 2 α m 2 N A 2 N A 1 - m 1 N A 1 α N A 2 -1 2 α m 2 N a 2 N a 1 - m 1 N a 1 α N a 2 N a 1 N A 1 +N a 1 - N a 2 N A 2 +N a 2 2 - α m 2 N A 2 N A 1 + m 1 N A 1 α N A 2 - N a 1 N A 1 +N a 1 + N a 2 N A 2 +N a 2 4 N a 1 N A 1 +N a 1 + N a 2 N A 2 +N a 2 4 N A 2 N A 2 +N a 2 - N A 1 N A 1 +N a 1 2 N A 1 N A 1 +N a 1 + N A 2 N A 2 +N a 2 4 - α m 2 N a 2 N a 1 + m 1 N a 1 α N a 2 - N A 1 N A 1 +N a 1 + N A 2 N A 2 +N a 2 4       

D.1.1 Resolution.

Following Remark 1, we recall that we assume that no major-effect allele has fixed. Here, we show that there is at most one instantaneous ecological equilibrium at each Z-level (for Z ∈] -1 -η A +η a 2 , 1 -

η A +η a 2
[), thanks to Proposition D.1 and Proposition D.2.

Proposition D.1. Suppose that no major-effect allele has fixed. Then, for Z ∈]-1-η A +η a 2 , 1-η A +η a 2 [, S 0 (Z) has exactly one solution (N a 1 , N A 1 , N a 2 , N A 2 ) ∈ (R * ) 4 , given by:

N a 1 = Y A N 1 -N 2 Y A -Y a , N a 2 = Y a Y A N 1 -N 2 Y A -Y a , N A 1 = N 2 -Y a N 1 Y A -Y a , N A 2 = Y A N 2 -Y a N 1 Y A -Y a ,
where the quantities (Y A , Y a , N 1 , N 2 ) are defined by (29):

                             Y A = g1 α m2 η A + η a + 2 (Z + 1) η A -η a 2    1 + m1 m2 4 g1 g2 η A -η a 2 2 1-η A +η a 2 +Z 2 + 1    , Y a = g1 α m2 η A + η a + 2 (Z + 1) η A -η a 2    1 + m1 m2 4 g1 g2 η A -η a 2 2 1-η A +η a 2 +Z 2 -1    , N 1 = 1 -g 1 (Z + 1 + η A ) 2 -m 1 + α m 2 Y A , N 2 = λ -g 2 (Z -1 + η A ) 2 -m 2 + m1 α Y A . ( 29
)
This solution (N a 1 , N A 1 , N a 2 , N A 2 ) defines viable numbers of each allele in each sub-populations if and only if:

[Y A N 1 > N 2 ] and [N 2 > Y a N 1 ]. ( 30
)
Proof. Let us introduce the following change of variables, valid under the assumption that no majoreffect allele has fixed:

N 1 := N A 1 + N a 1 , N 2 := N A 2 + N a 2 , Y A := N A 2 N A 1 , Y a := N a 2 N a 1 .
Then, under the assumptions made in Remark 1, the system (S 0 (Z)) is equivalent to:

               α m 2 Y a -m 1 + 1 -N 1 -g 1 (Z + η a + 1) 2 = 0, α m 2 Y A -m 1 + 1 -N 1 -g 1 (Z + η A + 1) 2 = 0, m1 α 1 Y a -m 2 + λ -N 2 -g 2 (Z + η a -1) 2 = 0, m1 α 1 Y A -m 2 + λ -N 2 -g 2 (Z + η A -1) 2 = 0. ( 31
)
This is equivalent to the following system:

         α m 2 (Y a -Y A ) + g 1 η A + η a + 2 (Z + 1) (η A -η a ) = 0, m1 α ( 1 Y a -1 Y A ) + g 2 η A + η a + 2 (Z -1) (η A -η a ) = 0, N 1 -1 -g 1 (Z + 1 + η A ) 2 -m 1 + α m 2 Y A = 0, N 2 -λ -g 2 (Z -1 + η A ) 2 -m 2 + m1 α Y A = 0.
As Z ̸ = 1 -η A +η a 2 , the closed subsystem on (Y A , Y a ) is, in turn, equivalent to:

Y A -Y a = A 1 (Z) := g1 α m2 η A + η a + 2 (Z + 1) (η A -η a ), -Y A Y a = A 0 (Z) := g1 m1 α 2 g2 m2 η A +η a +2 (Z+1) η A +η a +2 (Z-1) .

Y A and -Y a are the roots of the polynomial: P (X) = X 2 -A 1 (Z) X + A 0 (Z).

P has two real roots of opposite signs if and only if:

[A 0 (Z) < 0] ,
which is equivalent to:

-1 - η A + η a 2 < Z < 1 - η A + η a 2 .
Under the last condition on Z, A 1 (Z) is positive, A 0 (Z) is negative and we get:

           Y A = A1(Z) 2 1 -A0(Z) A 1 (Z) 2 2 + 1 , Y a = A1(Z) 2 1 -A0(Z) A 1 (Z) 2 2 -1 , (32)
which is equivalent to (29).

Inverting the initial change of variables leads to: Proof. Using the notation N 1 := N A 1 + N a 1 and N 2 := N A 2 + N a 2 , we compute thanks to a symbolic computation tool (Mathematica©):

N a 1 = Y A N 1 -N 2 Y A -Y a , N a 2 = Y a Y A N 1 -N 2 Y A -Y a , N A 1 = N 2 -Y a N 1 Y A -Y a , N A 2 = Y A N 2 -Y a N 1 Y A -Y a , hence (29)
det(J S L ) = - 1 4 m 1 α N 1 N 2 + α m 2 N 2 N 1 + 2 m 2 1 α 2 N a 1 N A 1 N a 2 N A 2 + 2α 2 m 2 2 N a 2 N A 2 N a 1 N A 1 + 2m 1 m 2 N A 1 2 N a 2 N a 1 N 1 N 2 + N a 1 2 N A 2 N A 1 N 1 N 2 + N A 2 2 N a 1 N a 2 N 1 N 2 + N a 2 2 N A 1 N A 2 N 1 N 2 + 2 N A 2 N a 1 N 1 N 2 + 2 N a 1 N a 2 N 1 N 2 < 0.

D.1.2 Stability.

Convergence toward a limit system locally in time in a slow-fast analysis relies essentially on a stability criterion of the fast equilibria which constitute the slow manifold [START_REF] Levin | Singular Perturbations of Non-linear Systems of Differential Equations and an Associated Boundary Layer Equation[END_REF][START_REF] Dekens | Evolutionary dynamics of complex traits in sexual populations in a heterogeneous environment: how normal?[END_REF]. In this subsection, we show that all fast equilibria found in Proposition D.1 and Proposition D.2 for a level

Z ∈] -1 -η A +η a ∂ α m 2 N a 2 -m 1 N a 1 + N a 1 1 -(N a 1 + N A 1 ) -g 1 (Z + η a + 1) 2 ∂N a 1 = 1 -(N a 1 + N A 1 ) -g 1 (Z + η a + 1) 2 -m 1 -N a 1 = - α m 2 N a 2 N a 1 -N a 1 ,
for (N a 1 , N A 1 , N a 2 , N A 2 ) solves S 0 (Z). The same holds for the other diagonal entries. 2. Let: χ J S 0 (X) = X 4 -tr (J S0 ) X 3 + b X 2 + c X + det J S0 , be the characteristic polynomial of J S0 . Let us verify the Routh-Hurwitz criterion: all the eigenvalues of J S0 are located in the left open half plane if and only if:

(i) det J S0 > 0, (ii) -tr (J S0 ) > 0, (iii) -tr (J S0 ) b -c > 0, (iv) (-tr (J S0 ) b -c) c -tr (J S0 ) 2 det J S0 > 0.

We have:

det J S0 = m 1 m 2 N a 1 N a 2 -N A 1 N A 2 2 1 N a 1 N A 2 + 1 N A 1 N a 2 > 0.
and:

tr

(J S0 ) = N 1 + N 2 + √ m 1 m 2 α √ m 2 N a 2 √ m 1 N a 1 + √ m 1 N a 1 α √ m 2 N a 2 + α √ m 2 N A 2 √ m 1 N A 1 + α √ m 2 N A 2 √ m 1 N A 1 > 0.
With the help of a symbolic computation tool (Mathematica©), we verify that the left hand side of the two last conditions are sums of positive terms, but are too long to be displayed here.

The Jacobian matrix of the linear system S L (Z) is exactly J S L and we also show that J S L satisfies the Routh-Hurwitz criterion: Proposition D.4. J S L has all its eigenvalues located in the left open half plane.

Proof. Let the following be the characteristic polynomial of J S L : χ J S L (X) = X 3 -tr(J S L )X 2 -1 2 tr(J 2 S L ) -tr(J S L ) 2 X -det(J S L ).

We show that J S L satisfies the Routh-Hurwitz criterion:

(i) -det(J S L ) > 0, (ii) -tr(J S L ) > 0, (iii) 1 2 tr(J 2 S L ) -tr(J S L ) 2 tr(J S L ) + det(J S L ) > 0. We havedet(J S L ) > 0 from the proof of Proposition D.2 and:

-tr(J S L ) = 1 + √ m 1 m 2 α √ m 2 N a 2 √ m 1 N a 1 + √ m 1 N a 1 α √ m 2 N a 2 + α √ m 2 N A 2 √ m 1 N A 1 + α √ m 2 N A 2 √ m 1 N A 1 > 1 + 4 √ m 1 m 2 .
We verify that the l.h.s. of the last condition is a sum of positive terms.

E Proof of Proposition 3.1.

Let us define the quantities:

             Y A, * = 2 g η m 1 + m 2 4 g 2 η 2 + 1 , Y a, * = 2 g η m 1 + m 2 4 g 2 η 2 -1 , N * 1 = 1 -g η 2 -g -m + 2 g η 1 + m 2 4 g 2 η 2 , N * 2 = N * 1 .
(34) Proposition D.1 states that the latter defines a solution to S 0 (Z) given that Z = 0:

(N a, * 1 , N a, * 2 , N A, * 1 , N A, * 2 ) = N * 1 Y A, * -1 Y A, * -Y a, * , N * 1 1 -Y a, * Y A, * -Y a, * , N * 1 1 -Y a, * Y A, * -Y a, * , N * 1 Y A, * -1 Y A, * -Y a, * .
Since Y A, * > 1 and Y A, * Y a, * = 1, this solution is viable under the condition: N 1 > 0, hence requiring :

1 + 4 g 2 η 2 + m 2 > g η 2 + g + m, which in turn is equivalent to (13). Proposition D.2 next states that S L (Z) has a unique solution (δ * , δ A, * , δ a, * ) for such allelic population sizes (N a, * 1 , N a, * 2 , N A, * 1 , N A, * 2 ). One can compute that:

δ A, * = δ a, * = g 1 + η + Y A, * (1 -η) m(1 + Y A, * ) , δ * = - 2g 1 + η -Y A, * 2 (1 -η) Y A, * .
Finally, one can verify that (N a, * 1 , N a, * 2 , N A, * 1 , N A, * 2 , δ * , δ A, * , δ a, * ) along with setting Z * = 0 is a solution of the last equation of (12).

  et al. (2005) for asexual populations and adapted recently to sexual populations with the standard infinitesimal model (Calvez, Garnier, and Patout 2019; Patout 2020; Garnier et al. 2022) and develop new convex analysis arguments. Throughout this section, the time dependency will be omitted for the sake of clarity. Presentation of the methodology. Almost two decades ago, Diekmann et al. (2005) introduced a methodology to determine the dynamics of the trait values around which trait distributions get concentrated as Dirac masses under the regime of small variance. This methodology has since been used successfully to study several evolutionary questions, initially for asexual models, where the diversity generated by mutations of small variance is modelled by a linear operator translating the distribution of mutational effects (Perthame and Barles 2008; Barles, Mirrahimi, and Perthame 2009; Mirrahimi 2017; Mirrahimi and Gandon 2020). It has recently been adapted to study sexually reproducing populations

Figure 1 :

 1 Figure 1: Illustration of the Hopf-Cole transform to study concentration phenomena. This transformation unfolds singular distributions n ε close to a Dirac mass (in purple), by defining more regular proxies: U ε (in green) such that n ε = 1 √ 2πε e -Uε ε 2. This figure suggests that, when ε vanishes, the limit U is regular and positive and cancels at the support of the limit measure n.

  al. 2022; Calvez, Garnier, and Patout 2019; Patout 2020; Dekens 2022) and define the proxies

Figure 2 :

 2 Figure 2: Illustration of the slow-fast variables Z ε ,δ ε , δ Aε and δ a ε (in red), introduced in (11). This figure displays a situation where the two major-effect alleles are segregating in both habitats in a symmetrical fashion. The graph represents the two local trait densities for each of the two alleles: ñA 1,ε , ñA 2,ε , ña 1,ε , ña 2,ε (the same color is for the same deme, and the same linestyle is for the same major-effect allele), as a function of the trait ζ = z + η A (resp. z + η a ), where z is the infinitesimal contribution and η A (resp η a ) is the effect of the major-effect allele. In red, we indicate a visualization of the new variables introduced in (11). Z ε is the mean infinitesimal part of the metapopulation, δ ε the spatial average of the local difference between the two allelic mean infinitesimal parts, δ A ε and δ a ε the spatial discrepancies in the mean infinitesimal parts per allele. Note the difference in notation between the trait densities ñA i,ε and the infinitesimal contribution densities n A i,ε (which are the ones used in the analysis), which are linked by n A i,ε (z) = ñA i,ε (z + η A ) (respectively z + η a for ña i,ε ).

Figure 3 :

 3 Figure3: Stability region of the symmetrical polymorphic equilibrium (in yellow), for four major effects η ∈ {0.5, 0.7, 1, 1.3} (recalling that θ = 1), when m (y-axis) and g (x-axis) vary in [0, 3]). This figure highlights the gain and loss of polymorphism with regard to increasing selection, which is not predicted by the one-locus model (abbreviated as OLM in the legend), according to which polymorphism is maintained below the extinction threshold represented by the dashed yellow line. The stable region (in yellow) becomes larger as η grows closer to 1, as the major-effect allele effects can then allow local adaptation to the two patches on their own, then shrinks again. The red crosses in Fig.3aindicate the parameters used for the individual-based simulations (see Section 3.2 and Fig.5and Fig.9).

Figure 4 :

 4 Figure 4: Phase lines of the limit equation (15), when the migration rate and the strong allelic effect are held constant (m = 0.1, η = 12 ) and the selection strength varies (the lighter the color, the stronger the selection). Solid curves indicate that the system is polymorphic, whereas dashed curves indicate that one major-effect allele has fixed. Every intersection of the horizontal black line and a solid colored curve with a negative (resp. positive) slope indicates a locally stable (resp. unstable) polymorphic equilibrium. The darker curve with weak selection g = 0.01 has a positive slope at Z = 0 (unstable), the following curves have a negative slope at Z = 0 (stable for selection between g = 0.18 and g = 0.86), and finally the lightest curves have a positive slope at Z = 0 (unstable for g ≥ 1.03), which is consistent with Fig.3a. Note that there exists additionally two mirrored asymmetrical polymorphic equilibria for g = 0.86, for Z ≈ ±0.5 (indicated by the red arrows), which were unsuspected prior to this numerical analysis.

  (a) Major-effect locus with polygenic background: weak selection (g = 0.1). (b) Control case without polygenic background: weak selection (g = 0.1). (c) Major-effect locus with polygenic background: intermediate selection (g = 0.5). (d) Control case without polygenic background: intermediate selection (g = 0.5).(e) Major-effect locus with polygenic background: strong selection (g = 1).(f) Control case without polygenic background: strong selection (g = 1).

Figure 5 :

 5 Figure 5: Variance at the major-effect locus across time for increasing selection (top to bottom: g = 0.1, 0.5, 1) at a fixed rate of migration (m = 0.8), with symmetrical initial subpopulation sizes. p denotes the local frequency of the major-effect allele A.The left panel is obtained with both a major-effect locus (η = 1/2) and a polygenic background of 200 loci, whereas only the major-effect locus is present in the right panel. For each subfigure, 20 replicates simulations were run per set of parameters, according to the setting explained in Section 3.2. In each subfigure, the solid line represents the median trajectory and the shaded area indicates the 0.2 and 0.8 quantiles. The dashed lines represent the median trajectories of the numerical resolutions of the deterministic model (1). This figure confirms that polymorphism of the major-effect locus is maintained only when selection is intermediate in strength (panel c) in presence of a polygenic background (left panel). The small figures embedded in each figure of the left panel represent the dynamics of the subpopulation sizes (N 1 and N 2 ). They highlight the qualitative difference between the loss of polymorphism at the major-effect locus with weak or strong selection, as fixation occurs without change in subpopulation size with weak selection (Fig.5a) and as subpopulations sizes become asymmetrical with strong selection (Fig.5e).

  Kirkpatrick 2001; Hendry, Day, and Taylor 2001; Dekens 2022 and for asexually reproducing populations, see Débarre, Ronce, and Gandon 2013; Mirrahimi 2017; Mirrahimi and Gandon 2020; Hamel, Lavigne, and Roques 2021).

Figure 6 :

 6 Figure 6: Summary of the complete analytical outcomes of the model, for η = 0.5 and varying migration (y-axis) and selection (x-axis). The figure combines the results obtained in Section 3.1 on the stability of the symmetrical polymorphic equilibrium with the results of the model of Dekens (2022) (equivalent to this model upon loss of polymorphism). For bounded selection, when migration increases, there is a threshold over which the polymorphism at the major-effect locus is lost due to the blending effect of migration (consistent with Yeaman and Whitlock 2011). The population then becomes equally maladapted to both habitats (generalist -symmetrical monomorphism, in the green region). For this specific major-effect allelic effect η = 0.5, there exists additionally an interval of selection strength (≈ [0.7, 1]) for which the major polymorphism might not be stable at all migration rates below the critical threshold. This phenomenon does not seem to hold when the major effect is larger (see Fig.3b). For bounded migration (below the threshold rate over which the strong migration blending hampers the major polymorphism), when selection strength increases, the polymorphism at the major-effect locus (yellow region) is lost, and the population becomes adapted to one of the two habitats (specialist -asymmetrical monomorphism, in the blue region). As this figure is obtained in the small segregational variance regime (which should remain smaller than the other parameters of the system for the analysis to be valid), we warn that the outcomes displayed in the vicinity of the x-axis (very small migration rates) might not correspond to the limit when the migration rate is 0 (no migration).

0)

  Scaling of timet := ε 2 T (ε 2 := σ 2 min|a (k) | 2 ≪ 1 ⇝ few diversity via inf. model of reproduction)1) Formal analysis (justify Gaussian distributions -Proposition B.1):

Figure 7 :

 7 Figure7: Layout of the slow-fast analysis in Appendix D. This figure presents the key elements of the separation of time scales leading from P ε to P 0 . The stability of the fast equilibria (studied in the two subsystems S 0 (Z) and S L (Z)) is the crucial argument underlying the convergence result stated in Theorem 2.1. The resolution of S 0 (Z) leads to a uniqueness result that is unexpected with regard to the analogous resolution in Dekens (2022) (see Remark 5).

  . It defines a viable solution with positive entries if and only if YA N 1 > N 2 and N 2 > Y a N 1 . Proposition D.2. For all allelic sizes of subpopulations (N a 1 , N A 1 , N a 2 , N A 2 ) ∈ R * + 4and Z ∈ R, there exists a unique solution (δ, δ A , δ a ) to the system S L (Z).

  (a) Major-effect locus with polygenic background: weak selection (g = 0.1). (b) Major-effect locus with polygenic background: intermediate selection (g = 0.5).(c) Major-effect locus with polygenic background: strong selection (g = 1).

Figure 8 :

 8 Figure 8: Same as the left panel of Fig. 5, but with fixed subpopulations sizes.

Figure 9 :

 9 Figure 9: Same as Fig. 5, but with asymmetrical initial subpopulation sizes.

Figure 10 :

 10 Figure 10: Same as Fig. 5, but with a quantitative background of L = 50 loci and σ LE = 0.2.

Figure 11 :

 11 Figure 11: Same as Fig. 5, but with a quantitative background of L = 50 loci and σ LE = 0.2 and asymmetrical initial subpopulations sizes.

Acknowledgements

The authors thank John Wakeley and two anonymous reviewers for very constructive and detailed comments that, in their opinion, greatly improved the manuscript. L.D also thanks Sepideh Mirrahimi for valuable discussions and comments, and Florence Débarre for insightful conversations. L.D has received partial funding from the ANR project DEEV ANR-20-CE40-0011-01 and a Mitacs Globalink Research Award. S.P.O. has received funding from the NSERC Discovery Grant: RGPIN-2016-03711. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programm (grant agreement No 865711).

F Supplementary IBS with fixed subpopulations sizes

In this appendix, we show in Fig. 8 the analogous results as those presented in Fig. 5, but with a slightly different procedure for the IBS, which adjusts the birth rates to compensate exactly for the deaths by selection at each generation, thus keeping the subpopulations sizes fixed. As mentioned in the main text, the loss of the polymorphism at the major-effect locus still occurs with weak selection, but not with strong selection.

G Supplementary IBS with asymmetrical initial conditions or

different parameters for the genetic architecture (L = 50 loci and σ LE = 0.2)

In this appendix, we first show that the phenomenon of loss of polymorphism in the presence of a quantitative background with weak or strong selection is robust to asymmetrical initial population sizes, and even occurs much faster (Fig. 9). We emphasize on the excellent agreement of the deterministic iterations with the individual-based simulations Furthermore, we also show that our findings hold when considering a smaller number of loci involved in the quantitative background (L = 50 instead of 200), with increased relative effect (σ LE = 0.2 instead of 0.1), so that the trait range

] extends beyond the local optima (-1 and 1) even in the absence of major effects. We display the results of the IBS with symmetrical initial subpopulation sizes in Fig. 10 and with asymmetrical initial subpopulation sizes in Fig. 11. Note that the right panel of each figure does not change from Fig. 5 and Fig. 9, because the control case does not depend on the number of loci, but we choose to display it anyway for consistency of comparison. One can notice that the time to fixation at the major-effect locus in the presence of a quantitative background under weak (Fig. 10e, Fig. 11a) and strong selection (Fig. 10e, Fig. 11e) is reduced compared to when the quantitative background comes from a larger number of loci (Fig. 5, Fig. 9). Moreover, the sensitivity of the numerical resolutions of (1) with regard to symmetrical initial states is more pronounced here.