Polygenic barriers to gene flow: the role of dominance, haploid selection and heterogeneous genetic architectures

We study the maintenance of polygenic local adaptation and its effects on reproductive isolation in a mainland-island model for populations with a general biphasic life cycle, encompassing haploid and diploid models as special cases. We quantify the strength of a multilocus barrier to gene flow due to divergent local adaptation at L unlinked and weakly selected loci, and obtain predictions for the equilibrium frequencies of locally adaptive alleles with arbitrary dominance and haploid-phase selection, accounting for genetic drift using a diffusion approximation. We extend classical single locus results on the role of dominance in the mainland-island model to the multilocus case, highlighting how linkage disequilibrium in multilocus barriers has rather different effects on differentiation and swamping thresholds for recessive alleles compared to dominant ones. Details about the biphasic life cycle can be captured through a set of effective parameters, and we show that for the same total strength of selection over the life cycle, increasing the relative intensity of selection in the haploid phase leads to stronger barriers to gene flow. We study the effect of heterogenous genetic architectures of local adaptation on the resulting barrier to gene flow, characterizing the realized genetic architecture at migration-selection balance for different distributions of fitness effects. Our results highlight the importance of barrier heterogeneity in shaping observable patterns of differentiation between populations under divergent selection pressures.


Introduction
When a population is subdivided across multiple habitats with different environmental conditions, the extent to which distinct subpopulations can maintain locally beneficial genetic variation depends on the rate of migration between them. Migration between populations that maintain divergently selected alleles can generate migration load (a reduction in mean fitness due to the influx of locally maladaptive genes) or may lead to loss of local adaptation altogether (so-called swamping by gene flow) (e.g. Lenormand 2002). While local adaptation may be driven by a few conspicuous loci (e.g. adaptive melanism in peppermoths (van't Hof et al. 2016) or pocket mice (Nachman et al. 2003)), it is believed to typically be polygenic, involving alleles of different effect at many loci across the genome (Pritchard and Di Rienzo 2010;Le Corre and Kremer 2012;Westram et al. 2018; Barghi et al. 2020;Bomblies and Peichel 2022;Stankowski et al. 2023). When local adaptation is polygenic, migration from a population adapted to different environmental conditions will generate linkage disequilibria (LD) among selected loci, and the rate at which each individual invading locally deleterious allele is eliminated will be affected by such associations, a phenomenon often referred to as a 'coupling' effect (Barton 1983;Kruuk et al. 1999;Feder et al. 2012;Yeaman 2015;Sachdeva 2022). Such coupling effects will in turn affect the equilibrium migration load and swamping thresholds for the loci under selection. Neutral variation may also come to be associated with locally selected alleles, so that the latter constitute a 'barrier' to neutral gene flow, increasing neutral genetic differentiation (as quantified by F ST for instance) beyond the single locus neutral expectation (Bengtsson 1985).
Barrier effects due to divergent local adaptation at many loci may play an important role in the evolution of reproductive isolation, and hence speciation (Nosil 2012;Barton 2020). The colonization of a new habitat will often involve selection on polygenic traits and give rise to a subpopulation that exhibits some divergence from its ancestors (Barton and Etheridge 2018). Conditional on the initial succesful establishment of such a divergent subpopulation, whether or not speciation ensues depends on whether local adaptation can be maintained in the face of maladaptive gene flow (if any), and on the extent to which the partial reproductive isolation deriving from local adaptation may promote further divergence and strengthen reproductive isolation, either through reinforcement, coupling with intrinsic incompatibilities, or the establishment of additional locally beneficial mutations (Barton and De Cara 2009;Bierne et al. 2011;Butlin and Smadja 2018;Kulmuni et al. 2020). In this paper, we focus on the conditions under which polygenic local adaptation can be maintained in the face of maladaptive gene flow, assuming local fitness to be determined by an additive trait under directional selection.
Despite mounting evidence that local adaptation is indeed often polygenic (Bomblies and Peichel 2022), little is known about the underlying genetic details: How many loci are involved? What are the typical effect sizes? Are locally beneficial alleles typically closely linked or spread all over the genome? How non-additive is local adaptation? etc. (e.g. Yeaman and Whitlock 2011;Yeaman 2015;Bomblies and Peichel 2022). Moreover, even if such details were known, a number of open questions remain as to how the genetic architecture of local adaptation affects the ability of a population to maintain reproductive isolation in the face of gene flow. For instance, it is not directly obvious whether a heterogeneous architecture comprising both loci of small and large effect would admit more adaptive differentiation than a homogeneous architecture with the same average selective effect per locus. The closely related question of how much scope there is to infer the detailed genetic architecture underlying local adaptation from observed patterns of genomic differentiation, as for instance obtained through so-called 'genome scans', also remains largely unanswered. So far, most theoretical models have assumed rather simple architectures, dealing with biallelic loci with additive and equal effects on fitness (ignoring dominance and epistasis), that are either unlinked or uniformly spread along a block of genome (Barton 1983;Fraïsse and Sachdeva 2021;Sachdeva 2022); and statistical approaches for the inference of gene flow across the genome either make similarly crude assumptions (Aeschbacher et al. 2017), or ignore the genetic details of local adaptation altogether (Roux et al. 2013;Laetsch et al. 2022).
In a recent paper, Sachdeva (2022) showed that, when the loci under selection are unlinked, the effects of LD on equilibrium differentiation at any individual locus in a multilocus barrier can be well described by classical (deterministic or stochastic) single locus population genetic theory, provided that the migration rate m is substituted by an effective migration rate m e (Petry 1983;Bengtsson 1985; Barton and Bengtsson 1986;Kobayashi et al. 2008), which captures how gene flow at a focal locus is affected by selection against the associated genetic background. The effective migration rate for a neutral locus can furthermore serve as a quantitative measure of reproductive isolation (RI), i.e. RI = 1 − m e /m (Westram et al. 2022). Crucially, m e depends itself on the frequencies of divergently selected alleles, giving rise to feedback effects where a small increase in migration rate may cause a sudden collapse of local adaptation (i.e. swamping). In her paper, Sachdeva (2022) conducted a detailed study of the joint effects of drift and LD on swamping thresholds and neutral differentiation in the mainland-island and infinite-island models of population subdivision, assuming a haploid sexual life cycle and divergently selected loci of equal effect. In this paper, we extend the theoretical framework outlined in Sachdeva (2022), deriving an expression for the effective migration rate for a polygenic genetic architecture with arbitrary fitness and dominance effects across loci (referred to as a heterogeneous architecture or barrier) in a population with a haplodiplontic life cycle (which includes haplontic and diplontic life cycles as special cases). We use this m e to build an approximation for the marginal allele frequency distributions at migration-selection balance in a mainland-island model, and use these tools to assess how the effects of multilocus LD (i.e. coupling effects) on swamping thresholds and the genetic architecture of reproductive isolation depend on dominance, drift, life cycle assumptions and heterogeneity in fitness effects across loci.

Haplodiplontic mainland-island model
Here we outline a mainland-island model for a sexual population which may be subject to selection in both the haploid and diploid phases. We think of this model as a caricature of a bryophyte, pteridophyte, algal or fungal population, but as we shall see below, the model encompasses both diplontic and haplontic life cycles as well. Throughout, we shall assume that sexes need not be distinguished. We assume a regular and synchronous alternation of generations, where an island population of N haploids (gametophytes) produces an effectively infinite pool of gametes from which 2N k gametes are sampled that unite randomly to form N k diploid individuals (sporophytes), k being the number of diploids per haploid individual. The diploid generation produces in turn an effectively infinite pool of haploid spores through meiosis, of which N are drawn to form the next haploid generation. In each generation, we assume M haploid individuals on the island are replaced by haploid individuals from a mainland population, where M is Poisson distributed with mean N m. Fitness on the island is determined by L unlinked biallelic loci which are under divergent selection relative to the mainland. The mainland population is assumed to have a constant, but arbitrary, genetic composition. Unless stated otherwise, we shall assume the mainland to be fixed, at each locus, for the locally deleterious allele on the island. Fitness effects are allowed to vary arbitrarily across loci. Denoting the alleles at locus i by A i,0 and A i,1 , we designate by w i,j the relative fitness of the haploid genotype A i,j and w i,jk the relative fitness of diploid genotype A i,j A i,k . We suppress the index i when considering a generic locus. We assume throughout that w 0 = 1 and w 1 = e s 1 for the haploid phase, and w 00 = 1, w 01 = w 10 = e s 01 , and w 11 = e s 11 for the diploid phase. Throughout, we denote the frequency of the allele with relative fitness 1 (on the island) at locus i by p i , and the frequency of the alternative allele by q i = 1 − p i . Fitness is determined multiplicatively across loci, so that, for instance, the log relative fitness of a haploid individual fixed for all the '1' alleles (genotype A 1,1 , A 2,1 , . . . , A L,1 ) is given by log w = L i=1 s i,1 . We assume that each haploid (diploid) individual contributes gametes (spores) to the gamete (spore) pool proportional to its fitness. We assume symmetric mutation at a small constant rate u per locus, occurring at meiosis.
Individual-based simulations of this model are implemented in a Julia package (Bezanson et al. 2017) available at https://github.com/arzwa/MultilocusIsland. In the following sections, we build up a theoretical approximation to this fairly general multilocus model, and validate the approximations by comparing numerical results against individual-based simulations. We first derive the dynamics at a single locus, considering both deterministic and stochastic models. Next, we derive an approximation to the effective migration rate for the multilocus model using a rather general argument based on the reproductive value of migrant individuals. Lastly, we approximate the allele frequency dynamics of the multilocus model by plugging in the effective migration rate, which captures the effect of LD among selected alleles on the dynamics at a neutral locus, in the single locus theory.

Deterministic dynamics
We first consider a deterministic model for the allele frequency dynamics at a single locus, ignoring the influence of the other loci as well as genetic drift. As shown in detail in sec. S2.1, for weak selection and migration, the dynamics of p can be described in continuous time by the nonlinear ordinary differential equation (ODE) where p * is the frequency on the mainland of the allele which is beneficial on the island, s a = s 1 + s 01 and s b = s 11 − 2s 01 , the latter being a measure of dominance (i.e. the deviation from multiplicative fitnesses, sometimes called ι (Otto 2003;Manna et al. 2011)). Usually, s 1 , s 01 and s 11 will be assumed to be negative, and p * will be assumed to be small, so that selection increases p, whereas migration decreases p. When s 1 = 0, we obtain the standard diploid mainland-island model, which is commonly parameterized in terms of a dominance coefficient h and selection coefficient s, so that s 01 = sh and s 11 = s. When s 1 ̸ = 0 (i.e. there is selection in the haploid phase), we can define an effective selection coefficient s e = 2s 1 + s 11 and dominance coefficient h e = s 1 +s 01 2s 1 +s 11 (except when s e = 0), so that, under weak selection, the single locus dynamics for an arbitrary haplodiplontic life cycle can be described by a diploid model with these effective parameters. The equilibria of eq. 1 are analyzed in detail in sec. S2.2.

Diffusion approximation to the stochastic dynamics
Still considering a single locus, we now account for the effects of drift. Denoting by X n and Y n the number of A 1 copies in the nth haploid, respectively diploid, generation, the life cycle as outlined in sec. 2.1 corresponds, for a neutral locus, to the following Markov chain model: Note that one unit of time corresponds to a single alternation of generations, involving two sampling stages: first we sample 2N k gametes from an infinite pool of gametes with allele frequency X n /N for the A 1 allele (eq. 2), next we sample N haploid genotypes from an infinite pool of haploid spores where the allele is at frequency Y n /2N k (eq. 3). This model is akin to the standard Wright-Fisher (WF) model with variable population size, regularly alternating between N and 2N k gene copies. The corresponding effective population size is hence N e = (N −1 + (2N k) −1 ) −1 , twice the harmonic mean of the phase-specific number of gene copies (Hein et al. 2004) (twice because our unit of time is an alternation of generations, not a single generation).
A similar Markov chain model can be defined for a selected locus with migration and mutation, but we will not consider this explicitly here, immediately considering a diffusion approximation instead. When evolutionary forces are sufficiently weak, diffusion theory can be applied to approximate the equilibrium allele frequency distribution implied by such a Markov chain by a continuous density ϕ(p) of the form Felsenstein 2005), where for the haplodiplontic mainland-island model, the infinitesimal mean and variance will be, respectively, where we assume mutations occur with rate u for both alleles. This yields a probability density function for the equilibrium allele frequency distribution where no closed form expression is known for the normalizing constant. This is essentially Wright's distribution, generalized to a haplodiplontic life cycle (Wright 1937).

Effective migration rate
We now derive an expression for the effective migration rate m e , which captures the reduction in gene flow at a focal locus embedded in a genotype with multiple selected loci due to LD.
As shown formally in Kobayashi et al. (2008), for weak migration, the reduction in gene flow relative to the 'raw' migration rate m, termed the gene flow factor (gff), depends on the expected reproductive value (RV) of migrants in the resident background (i.e. the expected long-term contribution of a migrant individual to the neutral gene pool on the island, relative to individuals in the resident island population). At any time, the proportion of individuals with recent migrant ancestry on the island is O(m), so that the probability of individuals with migrant backgrounds mating with each other to produce, for instance, F2 crosses of the migrant and resident genotypes, is O(m 2 ), and hence negligible for sufficiently weak migration. The descendants of a migrant individual will therefore most likely correspond to F1 and subsequent backcross generations, so that to a good approximation, the RV of a migrant individual depends on the relative fitnesses of F1, BC1, BC2, etc. individuals.
denote the relative fitness of an individual derived from an nth generation haploid, respectively diploid, backcross of an initial migrant individual with the resident population (i.e. W (1) d is the relative fitness of an F1 diploid, W (2) d of an offspring from a F1 × resident cross (BC1 generation), etc.). Assuming migration occurs in the haploid phase before selection, the gff can be expressed as where W h is the relative fitness of the haploid migrant in the resident population (Barton and Etheridge 2018;Sachdeva 2022;Westram et al. 2022). Note that this involves an expectation over all possible lines of descent of an initial migrant spore. In practice, g is determined only by the first 10 backcross generations or so, as subsequent backcrosses are essentially indistinguishable from residents. In order to derive a useful approximate expression for g, we shall make two further important assumptions: (1) both the resident and migrant gene pool, as well as each backcross generation, is in Hardy-Weinberg and linkage equilibrium (HWLE); (2) the expected allele frequency at any locus in any backcross generation is midway between that of the parents (e.g. the mean of the mainland and island allele frequencies for the F1 generation). In reality, due to Mendelian segregation, individuals inherit not exactly half of the selected alleles of each parent, and this segregation variance will lead to variation within F1s, BC1s, etc. on which selection can act. This will cause deviations from the midparent value which are O(s 2 ), so that this last assumption becomes more plausible when local adaptation is due to more and more loci of smaller effect.
Under these assumptions, each of the W (n) is determined solely by the frequencies of the selected alleles in the mainland and the island populations at the assumed equilibrium. This allows us to determine E[q (n) i ], the expected frequency of the locally deleterious allele (in the island) at locus i among nth generation descendants from a migrant, in terms of the allele frequencies in the mainland and island population. Indeed, assumption (2) implies the recursive relation q i.e. the average number of selected alleles carried by an nth generation backcross is the mean of the number of such alleles carried by an (n − 1)th generation backcross and a resident individual. Hence, we have E[q (n) i ] = 1 2 n (q * i + (2 n − 1)E[q i ]). Denoting the selection coefficient at locus i for the haploid phase by s i1 , the expected relative fitness of an nth generation haploid descendant is where we have assumed that per-locus selection is sufficiently weak that O(s 2 ) terms can be ignored. For the diploid phase, a similar argument shows that for the (n + 1)th generation, where s i01 and s i11 are the selection coefficients against heterozygotes and homozygotes at locus i respectively, and where, analogous to the single locus model, Putting everything together, the approximate gff becomes where, similarly, s i,a = s i1 + s i01 . It is worth stressing that the gff is a function of the differentiation between the mainland and island population as well as the heterozygosity E[pq] on the island, and that, although we assume migration is sufficiently rare, we do not assume that alleles introduced by migrants are rare. We shall often highlight the dependence of the gff on the allele frequencies and heterozygosities by writing g d ]) 2 , i.e. the product of the relative fitness of a haploid migrant in the haploid resident population and the relative fitness of the diploid F1 in the diploid resident population, squared. Hence, g can, in principle, be determined empirically.
If we assume all loci to have the same selection and dominance coefficient (a homogeneous barrier), and that the mainland is fixed for the locally deleterious allele on the island, eq. 6 can be simplified to where we have expressed s a = −s e h e and s b = −s e (1 − 2h e ) in terms of the effective selection coefficient s e against the invading allele, and the effective dominance coefficient h e of the invading allele over the locally beneficial one (see sec. 2.2.1). Here, the first factor is just the gff associated with a haploid L-locus system with selection coefficients s e h e . The second factor captures the effects of dominance and depends on the heterozygosity E[pq]. Clearly, h e has opposing effects on both factors. The immediate effect of dominance is therefore that the gff is decreased (barrier strength increased) relative to the additive case whenever invading alleles exhibit a dominant deleterious effect on the island (h e > 1/2). Only when heterozygosity (E[pq]) becomes appreciable does the second factor contribute to the increase (when h e > 1/2) or decrease (when h e < 1/2) of the gff. The implications of these observations for the maintenance of adaptive differentiation will be explored in detail in the results section.
Two remarks are due. Firstly, the gff as derived above yields the effective migration rate at an unlinked neutral locus. We can calculate the gff at a selected locus by making the assumption that it is the same as that of a hypothetical neutral locus at the same locationan assumption which is only expected to work well if selection at the focal locus is sufficiently weak. Hence, if we wish to calculate the effective migration rate for a selected locus in the barrier, say locus j, the relevant gff is obtained by excluding index j from the sum in eq. 6. Secondly, as in the model outlined in sec. 2.1, we have assumed that migration occurs at the start of the haploid phase, reflecting a process such as spore dispersal. However, it should be emphasized that, while the details of when migration occurs in the life cycle do not matter for the single locus model as long as selection and migration are sufficiently weak (so that the continuous-time limit is appropriate), these details do matter for the effective migration rate. This is because, although selection per locus is weak (s being small), selection against migrant genotypes can be strong (Ls being appreciable). Thus, when migration is due to dispersal of gametes (e.g. pollen dispersal), the first generation experiencing selection on the island will be the diploid F1 generation, so that the appropriate gff under the same approximation is Secondly, when migration occurs at the beginning of the diploid phase (e.g. seed dispersal), the first generation experiencing selection will consist of diploid migrant individuals, so that gE[W If the haploid, diploid and gametic migration rates are m 1 , m 2 and m 3 respectively, the effective migration rate will be (m 1 + E[W h ] −1 m 3 )g. Unless stated otherwise, in the present work, we shall assume migration to be due to dispersal of haploid spores, so that eq. 6 gives the relevant gff.

Dynamics and equilibria for the multilocus model
The gff captures the effect of LD among selected loci on the rate of gene flow from the mainland into the island at any individual locus. The key observation is that a certain separation of time scales applies: although selection against migrant genotypes can be very strong in the polygenic case (of magnitude Ls, roughly), selection at any individual locus is still assumed to be weak, so that, when linkage is weak or absent, LD among selected loci becomes negligible after an evolutionarily short period in which entire sets of alleles are efficiently removed together. Hence, on the longer time scales at which migration-selection balance is attained, the allele frequency at any individual locus should essentially follow the single locus dynamics, with LD reducing the effective migration rate by a factor equal to the gff (Sachdeva 2022). As a consequence, in the deterministic case, we expect that the effects of LD should be well captured by substituting the effective migration rate m e = mg for m in eq. 1. Specifically, we get a system of L coupled differential equations, where for 1 ≤ j ≤ L, where we assumed the mainland to be fixed for the deleterious allele on the island at all loci. Here we write g j [p −j ] for the gff as in eq. 7, to highlight the dependence of the gff at locus j on the allele frequencies at the other L − 1 loci. Note that in the deterministic model, the expected values in eq. 6 disappear, i.e. at any time E[q j ] = q j and E[p j q j ] = p j q j . We study the equilibria of this model by numerically solving for p at stationarity (ṗ j = 0, for 1 ≤ j ≤ L).
As in Sachdeva (2022), we can also plug m e into the single locus diffusion approximation to determine the equilibrium allele frequency distribution for each locus on the island. Specifically, we postulate that the joint distribution of allele frequencies factorizes as where Z is a normalizing constant and ϕ was defined in eq. 4. Eq. 9 can be thought of as the distribution associated with a Markov random field over the complete graph with L vertices. The marginal allele frequency distribution at any particular locus depends on the allele frequencies at the other L − 1 loci. We can compute moments of the allele frequency distribution at each locus by solving the whole system self-consistently, that is, by assuming where the Z's are again normalizing constants, and E[pq −j ] is the vector of expected heterozygosities at all loci excluding locus j ). To solve this system of 2L nonlinear equations, we use the fixed point iteration outlined in sec. S2.3. The numerical methods used in this paper are also implemented in the Julia package available at https://github.com/arzwa/MultilocusIsland.

Realized genetic architecture of local adaptation
To determine the realized genetic architecture of local adaptation at migration-selection balance for heterogeneous barriers, we calculate the conditional probability density for the selection and dominance coefficient at a locus, given that a divergent allele is observed on the island for that locus, i.e.
where f DFE denotes the joint density of the selection and dominance coefficient in the L-locus barrier, X i is an indicator random variable (equalling 1 when a randomly sampled allele on the island at locus i is of the locally beneficial type and zero otherwise), B is a shorthand for the selection and dominance coefficients at the L − 1 other loci ('B' for background), and we integrate over the set of all possible such backgrounds B. Note that f DFE is equivalent to f (s i , h i |X i = 1) in the absence of migration. For a given DFE model, we can characterize this conditional probability density using a Monte Carlo approach by sampling random L-locus genetic architectures from the DFE and calculating for each (s i , h i ) pair in the barrier the expected beneficial allele frequency E[p i |s i , h i , B] as a weight. The weighted sample will be distributed according to f .

Results
The results are organized as follows: we start with an analysis of the deterministic multilocus system, examining the effects of LD on equilibrium allele frequencies and contrasting the roles of dominance in the single locus model with the multilocus model, neglecting drift and heterogeneity of selective effects across the barrier. We next consider the effects of genetic drift and investigate how equilibrium allele frequencies on the island depend jointly on the population size, the strength of migration relative to selection at a single locus, the extent of LD, and dominance. We also consider the role of selection in both the haploid and diploid phase of a biphasic life cycle, showing how life cycle details can be accurately dealt with using a set of effective parameters. Lastly, we study the effects of heterogeneous barrier architectures on the maintenance of local adaptation and patterns of equilibrium differentiation.

Multilocus barriers with dominance in the deterministic model
We first analyze the deterministic multilocus model for a homogeneous barrier. We shall assume the mainland to be fixed for the locally deleterious allele (in the island habitat) at all loci. We only consider diploid selection here, with s 01 = −sh and s 11 = −s, where s is the selection coefficient against the locally deleterious allele, and h = s 01 /s 11 is the dominance coefficient. We emphasize that h measures dominance of the mainland (invading) allele over the island (locally beneficial) allele, so that h = 1 corresponds to a situation where the invading allele is fully dominant, or, equivalently, where the allele that confers local adaptation on the island is recessive. For the case where migration is at the haploid stage, restricting the analysis to diploid selection incurs no loss of generality, as haploid selection then simply amounts to a rescaling of the dominance and selection coefficients (see methods). The effective (diploid) dominance coefficient when there is haploid selection with intensity s 1 will be h e = s 1 +s 01 2s 1 +s 11 , so that the effect of haploid selection is to pull h e towards the additive case (h = 1/2) where selection acts on each gene copy independently, as it does in a strictly haploid model.

Effect of dominance on equilibrium frequencies
In the homogeneous deterministic model, all loci have the same dynamics if the initial allele frequencies are equal. From eq. 8, we find that in that case, the frequency p of the locally beneficial allele at any selected locus in an L + 1 locus system must satisfy at equilibrium We can solve this numerically for the equilibrium allele frequency. Note that if we set g[p] = 1, we recover the classical diploid single locus model, for which the equilibrium behavior is well understood (Haldane 1930;Nagylaki 1975). We briefly recapitulate the main results for the single locus model (see also sec. S2.2, and the gray lines in fig. 1). For the case h = 0.5 (no dominance, also referred to as codominance, or additivity), the equilibrium frequencyp of the locally beneficial allele decreases linearly from 1 to 0 as the rate of migration approaches the strength of selection per allele s/2. When local adaptation is due to a dominant allele (so that the invading allele acts recessively to reduce fitness on the island, i.e. h = 0), the migration rate beyond which no polymorphism can be maintained is increased to s, whilep is decreased relative to the additive case as long as m < s/4. When local adaptation is due to a recessive allele (h = 1), the model has two equilibria, one stable equilibriump + > 1/2 and one unstable equilibriump − < 1/2, as long as m does not exceed s/4. When this critical threshold is passed, swamping occurs for any initial frequency. Below this threshold, i.e. for m < s/4, whether or not polymorphism is maintained depends also on the history of the population: the island population cannot fix a new recessive beneficial variant, but an established recessive variant can be maintained upon secondary contact. Similar bistable behavior occurs for partial recessivity as long as h > 2/3.
We now consider equilibrium behavior in the multilocus case, where LD will cause deviations from single locus predictions. Fig. 1 shows the equilibrium behavior for a number of example parameter sets. As expected, stronger net selection against maladapted genotypes (stronger coupling, i.e. larger Ls) increases the equilibrium frequency of the locally beneficial allele relative to the single locus prediction, but the magnitude of this effect depends quite strongly on dominance. When invading alleles are recessive (h = 0; fig. 1 A), gene flow is not at all impeded when deleterious alleles are rare on the island (the gff being near one; fig. 1

E, F).
This is essentially because, irrespective of how many deleterious alleles a migrant carries, deleterious alleles will not be found in homozygotes as long as migration is sufficiently weak, and hence will not be 'seen' by selection. Only once deleterious alleles are segregating at appreciable frequencies on the island, are F1, BC1, etc. individuals likely to be homozygous at several loci, thus exposing (partially) recessive invading alleles to selection and reducing the RV of migrants. As a result, when invading alleles act recessively, a strong genome-wide barrier effect emerges only once differentiation falls below a critical threshold. The situation is clearly different when migrant alleles are dominant, as invading alleles will immediately express their full load in the resident population, irrespective of the alleles on the other haplotype, yielding efficient selection against migrant alleles (the gff being at its minimum when migrant alleles are rare, fig. 1 E, F). Any increase in the frequency of the deleterious allele on the island will merely increase the expected relative fitness of individuals with migrant ancestry, and hence reduce the efficiency of selection. We observe a transition between these two qualitatively different types of behaviour at intermediate values of h: when h < 1/3, the barrier strength, as measured by g −1 (Barton and Bengtsson 1986), increases as the locally deleterious allele increases in frequency on the island (and hence as differentiation between mainland and island decreases), decreasing the rate of gene flow, until a value of q = (3h − 1)/(4h − 2) is reached ( fig. 1, E, F).

Effect of dominance on swamping thresholds
In the single locus model, arbitrarily small frequencies of the locally beneficial allele can be maintained at migration-selection balance when h < 2/3, whereas in the case of h > 2/3, a sharp swamping threshold is observed as the locally beneficial allele reaches some critical frequency p c ≤ 1/2 (sec. S2.2, gray lines in fig. 1). Sachdeva (2022) showed that such sharp thresholds for swamping also appear in the absence of dominance due to coupling effects. LD both increases the critical migration rate (m c ) at which swamping occurs and the minimum level of differentiation that can be maintained before the swamping point is reached (p c ).
Our results indicate that dominance has a considerable influence on how LD sharpens and displaces swamping thresholds ( fig. 1 D, sec. S2.4). For h < 2/3 (i.e. when local adaptation is not strongly recessive), the critical migration threshold for swamping increases once Ls is sufficiently large (as in the additive case), albeit only marginally for moderate levels of divergence (Ls < 2, say) (sec. S2.4, fig. S2). This is in sharp contrast with the case where local adaptation is due to strongly recessive alleles (h > 2/3), where the threshold for swamping increases rapidly with Ls ( fig. 1 D). Importantly, the critical differentiation (p c ) level below which local adaptation collapses is strongly affected by dominance. In the additive case, one can show that sharp thresholds for swamping emerge as soon as Ls > 1 (sec. S2.4), in which case p c = 1 − 1/Ls, and hence arbitrary differentiation can be maintained near the critical point depending on Ls. For completely dominant local adaptation (h = 0), however, p c increases from 0 to a maximum of 1/2 as Ls → ∞, whereas for recessive local adaptation (h = 1), p c increases from 1/2 to 1 as Ls grows. This means, in particular, that for moderate levels of divergence, Ls > 0.75 say, and large population sizes, one would not expect to see locally beneficial recessives at frequencies much below 0.8, compared to 0.5 for the single locus model ( fig. 1 C). Fig. 1 further highlights the nontrivial feedbacks between the observed differentiation and dominance: whereas a consideration of the gff in a regime where migrant alleles are rare would suggest that swamping thresholds and equilibria depend roughly on

Accounting for drift
While the deterministic analysis points towards important effects of dominance on equilibrium differentiation and thresholds for swamping when local adaptation is polygenic, it is important to assess to what extent these carry over to finite populations. Indeed, Sachdeva (2022) showed that, for small values of N e s, the sharp thresholds for swamping predicted by deterministic multilocus theory need not apply, and that the critical migration rate may be significantly reduced. Furthermore, as long as selection is not much stronger than drift (i.e. N e s is not very large), the observed frequency of a locally adaptive allele (and hence differentiation between the mainland and island) at migration-selection-drift balance may deviate substantially from the deterministic limit, so that it becomes important to understand the distribution of allele frequencies at equilibrium on the island.  S4). Not only can we reliably obtain the expected frequency of alleles on the island, we also obtain very good predictions for the entire allele frequency distribution ( fig. 2 A). The sharp swamping thresholds observed in the deterministic model, in particular with recessive local adaptation (h = 1), correspond to strongly bimodal allele frequency distributions in the stochastic model. As described in more detail in sec. S2.5, this may render our numerical approaches sensitive to the assumed initial state of the island population. This sensitivity is itself biologically relevant, corresponding to assumptions on the (recent) history of the populations considered at equilibrium. Throughout, we shall assume a scenario of secondary contact, so that the island starts as fixed for the locally beneficial allele at all loci.
For a homogeneous genetic architecture, the swamping threshold at any individual locus depends on the total barrier strength (coupling strength) Ls, the dominance coefficient h, and the strength of selection per locus relative to drift N e s. Concomitantly, for any given m/s, these three key parameters will determine jointly whether significant adaptive differentiation is to be expected at a particular locus ( fig. 2 B). Unsurprisingly, the general consequence of genetic drift is to reduce the barrier effect, hence decreasing the expected differentiation at equilibrium ( fig. 2 B; fig. S5). Swamping thresholds are both decreased and made less sharp by drift, and the detailed behavior depends on the dominance coefficient. In the case with recessive local adaptation (h = 1), we see that sharp thresholds for swamping appear even when drift is quite strong (N e s > 4, say), and that, as Ls increases, the critical m/s value increases rapidly with increasing N e s ( fig. 2 B, bottom row). As expected from the deterministic theory, sharp thresholds for swamping in the additive and dominant case appear once Ls > 1, but only when drift is not too strong ( fig. 2 B, middle and top row). Considering, for instance, the predictions for Ls = 2 ( fig. 2 B, rightmost column), we see that for the additive case, a sharp threshold appears roughly when N e s > 10, whereas for the dominant case this only happens near N e s > 20. In the case where local adaptation is dominant, we now see clearly that when N e s is substantial (e.g. N e s = 16 in fig. S5), two phases can be distinguished as m/s increases. As noted above, the barrier strength initially increases with migration as deleterious alleles increase in frequency on the island, exposing more invading recessives to selection. However, above a certain critical deleterious allele frequency (E[q] > q c > 0.5), the barrier strength starts to decline again with increasing m ( fig. S5, h = 0 and h = 0.25 panels). This contrasts with the recessive case, where a positive feedback between the gff and the frequency of the deleterious allele generates a sharp threshold for swamping once a critical frequency q c < 0.5 is surpassed.

Selection in both the haploid and diploid phase
We now consider in more detail the case when selection acts both in the haploid and diploid phase. We parameterize the general model in such a way that we can investigate, for a given total barrier strength, the effects of the relative strength of selection in the diploid and the haploid phase and the degree of dominance (h) in the diploid phase. To this end, we assume where 0 ≤ τ ≤ 1 measures the relative strength of selection in the diploid phase (if one assumes selection to operate with constant intensity throughout the life cycle, this can be interpreted as the relative length of the diploid phase). Similar models have appeared in the study of life cycle modifiers (e.g. Otto 1994;Scott and Rescan 2017). Recall furthermore that we assume a regular alternation of N haploid and N k diploid individuals. Fig. 3 shows that we can accurately predict equilibrium allele frequencies for haplodiplontic life cycles with selection in both phases using the diffusion approximation, which depends on N, k, s, h and τ only through N e s e = N e (2 − τ )s and showing that, at least for weak selection, life cycle details can be accounted for by means of suitable effective parameters ( fig. S8). Fig. 3 suggests furthermore that, for a given total strength of selection Ls, predominantly haploid populations should be able to maintain more adaptive variation, and exhibit stronger reproductive isolation, irrespective of the degree of dominance in the diploid phase. This is because the effective selection coefficient in the model defined by eq. 12 is (2 − τ )s, so that the strength of selection per gene copy in haploids is twice that in diploids in the absence of dominance. Therefore, although we have shown in our analyses for diploids above that recessive local adaptation (h = 1) can lead to a stronger barrier to gene flow (if Ls is sufficiently large), increasing the relative strength of diploid selection (τ ) in the haplodiplontic model when locally beneficial alleles act recessively (thereby increasing h e ) does not strengthen the barrier, because this decreases s e at the same time. The relevance of these observations for the evolution and maintenance of haplodiplontic life cycles is however not very clear, as a life cycle modifier need not keep the overall strength of selection constant (Scott and Rescan 2017). Lastly, we note that the phase of the life cycle in which migration occurs can have an impact on the ability to maintain adaptive differentiation. When there is selection in both phases, and migration occurs in the diploid phase, maladapted diploid migrants are directly exposed to selection on the island, whereas this is not the case when migration occurs in the haploid phase. In the latter case, the first diploid generation exposed to selection is an F1 cross between the mainland and island, so that there is no selection on the homozygous diploid effect. We would therefore expect that migration in the diploid phase generally leads to stronger barriers to gene flow when selection acts in both phases, and this is indeed observed (fig. S9).

Heterogeneous genetic architectures
We now depart from the unrealistic assumption of equal effects and dominance coefficients across the polygenic barrier. In sec. 2.3, we developed the multilocus theory for potentially heterogenous (unlinked) genetic architectures, where the selection coefficients s 1 , s 01 and s 11 can vary arbitrarily across loci, and we verify that we do indeed obtain accurate predictions also in this setting ( fig. S10). This allows us to address in more detail a number of questions pertaining to the genetic architecture of local adaptation at migration-selection balance, accounting for both LD and genetic drift. Recall that the scenario we are considering is the following: at some point in time, the island is colonized by individuals from the mainland and rapid adaptation to the local conditions from standing genetic variation ensues (by driving L locally beneficial alleles to high frequencies on the island). After this initial idealized phase of directional selection on standing variation, we assume a situation of secondary contact, where, on average, N m haploid individuals on the island are replaced by mainland individuals in each generation, and the island evolves to an equilibrium state. We then ask what sort of loci can contribute to local adaptation and how the resulting barrier to gene flow leads to observable allelic differentiation. We also consider to what extent selective interference among loci results in departures from single locus predictions. We emphasize that we do not explicitly consider the buildup of divergence between the mainland and island population, nor some plausible model for what sort of standing variation is the source of the initial polygenic response, but merely ask how a given genetic architecture of local adaptation affects observable patterns of adaptive differentiation when there is gene flow. All results discussed in the remaining sections are obtained using the numerical approximations based on diffusion theory.

Effect of variation in fitness effects across a polygenic barrier
We first consider the case with variable selection coefficients across the L loci in the barrier, assuming no dominance. Fig. 4 (A) shows the average per-locus expected differentiation across the barrier (∆) when selection coefficients are sampled from a Gamma distribution (L = 100, Ls = 1). When migration is weak relative to selection (roughly m/s < 1/4), increasing the variance in fitness effects, while keepings constant, yields on average lower equilibrium differentiation than a homogeneous barrier of strength Ls, although the barrier strength (as measured by the average gff across the L selected loci,ḡ) is hardly affected ( fig. S11). At higher migration rates, where loci with selection coefficients close tos become prone to swamping, heterogeneous architectures tend to yield higher equilibrium differentiation and a stronger barrier effect than a homogeneous one with the same total effect Ls (figs. 4, S11). One should be careful, however, in the interpretation of∆. As shown in fig. 4 (B, C), differentiation across loci in the barrier often shows a strongly sigmoidal pattern, especially when Var[s] is large, where most loci are either strongly differentiated or not at all, and with rather few loci having E[p] near∆. This implies that empirically, instead of detecting L selected loci with an average differentiation of∆, we are more likely to observe about L∆ strongly differentiated loci. For low m/s, the weaker differentiation observed for more heterogeneous barriers is due to a smaller number of loci effectively contributing to local adaptation, with about half of the locally beneficial alleles swamped at m/s = 0.1 while the other half is strongly differentiated ( fig. 4 B), whereas the increased differentiation relative to the homogeneous case for larger m/s appears to be due to the presence of a larger subset of loci that resist swamping ( fig. 4 C). These results are not significantly affected when dominance coefficients are sampled independently from a symmetric Beta distribution ( fig. S12). The effect of increasing Var[h], while keeping the s i fixed across the barrier, is less dramatic than the effect of heterogeneity in selection coefficients, although we do see systematic increases and decreases in equilibrium differentiation depending on whether the migration rate exceeds the swamping threshold for recessives (which are associated with higher equilibrium frequencies) or not ( fig. S15).
Focusing on a single locus embedded within a heterogeneous barrier, we find that, for weak migration, variation in selection coefficients across the barrier has a negligible effect on differentiation at a focal locus with fixed selective effect, whereas (as already shown in fig. 4) it does have a strong effect on average differentiation across the L loci (figs. S13, S14). On the other hand, when migration is strong, a locus with selection coefficient s shows on average higher equilibrium differentiation when embedded in a heterogeneous barrier than in a homogeneous one, even when the average differentiation across the barrier is lower in the former (figs. S13, S14). These results are in line with the observed effect of Var[s] on g ( fig. S11) and indicate that the presence of a small number of loci of large effect in more heterogeneous barriers can have a strong effect on the expected differentiation at a focal selected locus ( fig. S13). As the distribution of selection coefficients is generally believed to be at least somewhat leptokurtic (note that excess kurtosis ∝ κ −1 for our Gamma DFE model), these results suggest that heterogeneity in selection coefficients can have important consequences for observable differentiation at migration selection-balance that would not be adequately captured by substituting an average selection coefficient in either single locus or multilocus theory.
Next we consider how barrier heterogeneity and selective interference affects different loci in the barrier differently, assuming selection and dominance coefficients to be independently distributed according to a Gamma and Uniform distribution respectively (see sec. S2.6). Comparing the expected differentiation at each locus in a heterogeneous polygenic barrier to the corresponding single locus predictions for that locus, we find, as expected, that the extent of selective interference is strongly dependent on the strength of migration relative tos and the total strength of selection Ls ( fig. 5), as well as the strength of genetic drift ( fig. S16). Differentiation is most strongly affected by the multilocus barrier effect for loci with recessive locally beneficial alleles, at least when m/s does not exceed the swamping threshold for those loci, whereas the deviation from the single locus prediction for dominant variants is considerably less (figs. 5, S13). Again we find that when migration is strong, increased heterogeneity of selection coefficients in the barrier generally leads to stronger selective interference, where an appreciable proportion of alleles are protected from swamping due to a few strongly selected barrier loci ( fig. S17; fig. S18). Fig. 5 further highlights that, even in the presence of a rather strong genome-wide barrier effect (e.g. Ls = 1.5), we would still expect to see considerable variation in equilibrium frequencies for locally selected alleles depending on their individual fitness effects.

The realized architecture of local adaptation at migration-selection balance
Although the above results indicate that, when Ls is appreciable, the effect of selective interference is strongest for recessive alleles ( fig. 5), this does not, however, imply that recessives necessarily contribute more to local adaptation at migration-selection balance than dominant alleles. Although strongly selected recessives will be associated with strong differentiation, weakly selected recessive alleles will be more prone to swamping than partially dominant ones (e.g. fig. S13). One way to quantify how these two phenomena interact to yield the realized genetic architecture of local adaptation (related to the concept of adaptive architecture, as defined in Barghi et al. (2020)) is by considering the conditional probability density for the selection and dominance coefficient at a locus, given that a divergent allele is observed on the island for that locus (see Methods, eq. 10). Fig. 5 (B) shows approximations to the marginal distributions f (s i |X i = 1) and f (h i |X i = 1) obtained in this way for the Figure 5: (A) Deviation of predicted expected allele frequencies for loci in heterogeneous polygenic barriers when accounting for LD (y-axis) from predictions based on single locus diffusion theory (x-axis). The rows show results for different total strengths of selection (different number of loci Ls withs = 0.01), whereas the columns show results for increasing rates of migration relative to selection (m/s). We assume the s i to be exponentially distributed with means and dominance coefficients are sampled uniformly from the [0, 1] interval. Each dot is associated with a single locus in an L-locus barrier, and is colored according to its dominance coefficient (yellow for locally beneficial recessives (h = 1), purple for dominants (h = 0)). Each plot shows results for 1000 such loci, subsampled from a total of 150000/L simulations of L-locus barriers. (B) Monte Carlo approximation to the marginal distribution of the selection and dominance coefficient conditional on observing a divergent allele on the island (i.e. f (s i |X i = 1) and f (h i |X i = 1), see eq. 10). The distribution graphed in gray shows f DFE , i.e. the marginal distribution of the selection and dominance coefficient for a random locus in the L-locus barrier in the absence of migration. We assumed N es = 20 and u/s = 0.005 for all results.
heterogeneous barrier model assumed in the preceding section.
As expected, we find that as migration rates go up (colors in fig. 5 B), the distribution of selection coefficients in the barrier at migration-selection balance shifts towards higher values of s, and that this effect becomes weaker with increasing Ls, which increases the extent by which small-effect alleles are protected from swamping. Notably, recessives contribute less to adaptation than dominants when migration is sufficiently strong, despite the fact that, conditional on no swamping, equilibrium frequencies of recessives are most affected by LD ( fig. 5 A). This is most notable when Ls is not large (top row in fig. 5). When Ls = 1.5 for instance, the depression in the conditional density at h = 1 becomes very slight even for relatively large migration rates (bottom row in fig. 5). We observe a similar shift in the distribution of dominance coefficients when h is Beta distributed with mean 2/3 instead of uniformly on the unit interval ( fig. S19). It is noteworthy that, despite s and h being Figure 6: The distribution of fitness effects (DFE) affects the realized genetic architecture of local adaptation at migration-selection balance. Contour plots for the joint density of h and s conditional on observing a divergent allele on the island (see eq. 10) are shown for the three DFE models (rows) for increasing rates of migration (columns). Values of∆ in the lower right corner denote the mean expected differentiation per locus. We assume Ls = 0.8,s = 0.01, N es = 20, u/s = 0.005 and exponentially distributed selection coefficients, and parameterize the DFE models so that E[h] = 2/3, assuming α = 2, β = 1 for the independent model, a = 7.2, b = 1.2, σ = 1 for the logistic model and K = 50 for the CK94 model (see sec. S2.6 for details on the different DFE models considered here). The densities are approximated using a Monte Carlo approach, simulating 500 replicate L locus genetic architectures from the assumed DFE model, determining the equilibrium allele frequencies for each replicate, and fitting a kernel density estimate to the sample so obtained. independent at each locus in the barrier, migration-selection balance induces a correlation between s and h in the distribution conditional on observed divergence, with variants of relatively large effect observed at equilibrium being more likely to act recessively than variants of small effect ( fig. S22, see also fig. 6, top row). The correlation is negligible for small migration rates, but as the strength of migration increases so that swamping effects become relevant, the correlation coefficient can become as large as 0.25, depending on Ls.
The simple DFE model assumed above where selection and dominance coefficients are independent is almost certainly inadequate. Theoretical and empirical work has indicated that, on the one hand, a correlation between selective effect and degree of dominance can be expected in the standing genetic variation that forms the basis for a polygenic selection response, with large-effect alleles more likely to act recessively (Caballero and Keightley 1994;Zhang et al. 2004;Agrawal and Whitlock 2011). On the other hand, it is well appreciated that during the process of adaptation, different loci enjoy different probabilities of rising to high frequencies, with dominant beneficial alleles having higher establishment probabilities than (partially) recessive ones with the same homozygous effect (Haldane's sieve;Haldane (1927), Turner (1981)). These two aspects interact when adaptation is from standing variation, as the on average higher initial frequency of partially recessive alleles increases the fixation probability, whereas its recessivity decreases it (Orr and Betancourt 2001). To examine how the realized genetic architecture of local adaptation depends on such assumptions, we consider two alternative, admittedly ad hoc, DFE models, outlined in sec. S2.6. Both models assume Gamma distributed selection coefficients and incorporate a positive correlation between s and h, so that alleles of large effect tend to be more recessive (recall once more that h in our case is the dominance coefficient of the invading allele, so h = 1 corresponds to recessive local adaptation). We keep the average dominance coefficient fixed to 2/3 for each model. In contrast with the independent model, we find that for the models that incorporate such a correlation between s and h, recessives are typically more likely to contribute to the realized differentiation at equilibrium (figs. 6, S19, S20, S21, S23, S24, S25). When we make the opposite assumption that locally beneficial alleles tend to be dominant (corresponding, for instance, to the case with a strong Haldane's sieve effect during adaptation), we find a somewhat less dramatic shift in the joint density as migration rates go up (fig. S26). The distribution also shifts towards higher selection coefficients, but somewhat less so than in the model with the opposite correlation. Swamping of partially recessive alleles of small effect further shifts the distribution towards smaller values of h. Clearly, these examples show how correlations between s and h among the loci under divergent selection can have a rather important influence on the realized genetic architecture at migration-selection balance (i.e. on which loci actually contribute to adaptive differentiation), driving up the relative contribution of recessives in one case but not in the other.

LD and polygenic migration-selection balance
Speciation, in essence, amounts to the buildup and maintenance of linkage disequilibria (Felsenstein 1981), with different sets of alleles maintained across populations. Heterogeneous selection can maintain different locally beneficial alleles across populations at multiple loci. When migration occurs between such divergently selected populations, sets of maladaptive alleles are introduced jointly, leading to statistical associations (i.e. LD) between selected alleles, and possibly strong selection against migrant alleles, in turn causing (partial) reproductive isolation (RI). These associations are broken down by recombination, so that the extent to which RI can be maintained depends on the relative strength of selection, migration and recombination. In this paper, we have focused on RI through the maintenance of polygenic local adaptation in the face of gene flow, assuming a scenario of secondary contact between two divergently adapted populations that have not diverged too much. The effects of polygenic selection against maladapted genotypes can be quantified by the gene flow factor g, or the associated effective migration rate m e = mg, which accounts for the (potentially quite strong) selection against individuals with recent migrant ancestry and the rapid decay of LD among surviving invading alleles (Barton and Bengtsson 1986). It should be noted that g has a phenotypic interpretation, and can, at least in principle, be experimentally determined.
We derived an expression for m e in a mainland-island model with heterogeneous barriers and a haplodiplontic life cycle, and showed how it can be used together with classical single locus population genetic theory to yield accurate predictions of equilibrium allele frequencies on the island. Importantly, this allows us to study the effects of coupling among barrier loci without assuming that locally deleterious alleles are somehow rare, enabling us to study swamping by gene flow in the polygenic setting. Our results show how the maintenance of adaptive differentiation in the face of gene flow depends jointly on the extent of LD, drift, dominance and variation in selective effects across loci. The general success of the approach indicates two important features of polygenic migration-selection balance. Firstly, it suggests that the 'separation of time scales' argument that is at the root of the approach indeed works, and does so beyond the haploid case with homogeneous selective effects (as treated in Sachdeva 2022). At least in the unlinked (and probably also weakly linked) case, strong selection against multilocus genotypes occurs only in the first couple of generations after a migrant arrives, and the long term fate of a migrant allele is unaffected by LD conditional on having survived these initial generations. As a consequence, the effects of LD are well described by the usual single locus dynamics, but with a reduced migration rate. Secondly, it indicates that our rather crude approximation to the expected reproductive value of a migrant individual on the island (which assumes HWLE within the island population, that migrants only cross with residents, and that in each such cross the proportion of migrant alleles is exactly halved) is an adequate estimator of the gff under weak migration. The approach enables us to study polygenic migration-selection balance in the mainland-island model using efficient numerical methods, and in particular to examine the relationship between the genetic architecture of locally adaptive traits and the ability to maintain adaptive differentiation in the face of gene flow.

Effects of dominance and haploid selection in polygenic barriers
Our analyses for homogeneous genetic architectures indicate that, when there is selection in the diploid phase, dominance can have a considerable impact on both the extent of adaptive differentiation and the ability to maintain it. Single locus theory can be somewhat misleading in this regard, in that it tends to emphasize the relative precariousness of local adaptation due to recessive alleles in the face of gene flow (because of the lower thresholds for swamping, e.g. Haldane 1930;Felsenstein 2005), whereas in the multilocus setting, depending on the total strength of divergent selection (Ls), partially recessive variants may lead to strongly increased swamping thresholds and produce a much stronger barrier to gene flow than dominant variants with the same homozygous effect ( fig. 1). The reason for this is that for recessive local adaptation, the dominant invading alleles are immediately exposed to selection, whereas for dominant local adaptation, the recessive invading alleles can introgress easily as long as the frequency of the locally deleterious allele is low. In the extreme case of maximal differentiation and completely dominant local adaptation, gene flow will be unimpeded whatever the extent of LD when migration occurs at the haploid stage. Our results show that the feedbacks between the level of differentiation and the strength of selection against migrants (as measured by the gff, which depends on the extent of differentiation) are strongly affected by the degree of dominance. It should be emphasized, however, that all our results assume a mainland-island model of migration and a scenario of secondary contact. The effects of dominance may turn out to be more subtle in models of population subdivision with multiple demes and more symmetric patterns of migration, in which case assumptions on environmental dependence of dominance may become important (e.g. Bürger 2013).
We developed our theory for a fairly general haplodiplontic life cycle, and showed that a model with selection in both the haploid and diploid phase can be brought into correspondence with a strictly diploid one through a set of effective parameters. Increasing the relative strength of selection in the haploid phase, for a fixed total strength of selection (Ls), yields stronger effective selection (larger s e ) at a selected locus, while making the fitness effect more additive (i.e. h e moves towards 1/2). One might therefore expect that local adaptation would lead to stronger barriers to gene flow in predominantly haploid species, but it is hard to make more refined predictions at such a general level. This relates to recent work on the strength of barriers to gene flow on sex chromosomes (Fraïsse and Sachdeva 2021) and in arrhenotokous species (Bendall et al. 2022), where increased exposure to haploid selection generally leads to stronger barrier effects. Importantly, similar considerations apply within the genome of haplodiplontic (and even diplontic) species, as there can be considerable variation in the relative expression levels in both phases across genes within a genome (e.g. Szövényi et al. 2011;Cervantes et al. 2023), so that it is likely that the relative strength of selection in the phases varies across the genome (Immler and Otto 2018). It therefore seems plausible that genes whose expression is biased towards the haploid phase may contribute more to local adaptation, similar to what has been described for sex chromosomes versus autosomes (Lasne et al. 2017), although this would of course depend on the relative extent of divergent selection in both phases. While general theoretical predictions are challenging to make without numerous additional assumptions, the fact that life cycle details can be accounted for using a few effective parameters can be useful in the empirical investigation of questions related to the relative importance of genetic variation in haploid-biased versus diploid-biased genes for local adaptation and reproductive isolation.

Heterogeneous architectures of polygenic barriers
When migration is not strong relative to the average strength of selection per locus, increased variation of s in the DFE underlying locally adaptive traits gives rise to lower overall differentiation at migration-selection balance, while generating a barrier to gene flow of similar strength (figs. 4, S11). On the other hand, for high rates of migration, a more heterogeneous architecture will tend to generate a stronger barrier to gene flow due to the presence of a larger number of strongly selected loci that resist swamping. As a consequence, increasing the variation among selection coefficients in a multilocus barrier (keeping Ls fixed) leads to stronger differentiation on average at a focal locus, especially at higher migration rates (figs. S13, S14). Hence, in a polygenic context, observed differentiation at any individual locus does not only reflect the fitness effect of that locus (i.e. s i and h i ), but may to a large extent be determined by a genome-wide barrier effect due to the presence of other divergently selected loci, depending on the number of barrier loci and their DFE, as well as the migration rate and population size. Nevertheless, our results also emphasize that, despite appreciable genome-wide coupling effects, considerable variation in equilibrium differentiation across non-swamped loci remains, depending on the allele-specific fitness effects ( fig. 5). Recent findings such as those of Stankowski et al. (2023) emphasize that readily discovered large-effect loci which bear a clear relationship to locally adaptive traits may often be associated with a polygenic background whose causal connection with local adaptation is less obvious. In line with this, our work suggests that when adaptive divergence (Ls) is appreciable, many loci under weak divergent selection that would be prone to swamping in isolation may persist due to the presence of a strong enough barrier.
Our findings concerning heterogeneous genetic architectures bear relevance to the inference of barriers to gene flow from genomic data. Recent approaches to quantify gene flow in pairs of diverging populations based on neutral variation have accounted for heterogeneity in barrier effects by assuming a demographic model in which the migration rate m (interpreted as m e ) varies along the genome. Roux et al. (2013) (see also Roux et al. 2014;) assume a model in which each locus is associated with a locus-specific m e distributed according to some genome-wide parametric family, and they infer the shape of this distribution using their ABC approach, whereas Laetsch et al. (2022) conduct (approximate) maximum likelihood inference of m e in windows along the genome, providing a model-based alternative to the popular F ST -based genome scans. Both approaches do not explicitly model aspects of the underlying genetic architecture of divergent selection that is assumed to cause variation in m e (and hence introgression probability) across the genome. So far, the only approach which has attempted to do so is the one proposed by Aeschbacher et al. (2017), where the authors combined information about local recombination rates together with deterministic population genetic theory (Aeschbacher and Bürger 2014) to obtain predictions of m e across the genome. They assume a set of unobserved selected loci with a fixed selection coefficient and no dominance to occur with a constant density across the genome, and infer the selection density per unit of map length (together with the migration rate) through its effect on observable neutral differentiation. We wonder, naturally, whether a similar inferential approach, but one allowing for drift and heterogeneous barrier architectures in the calculation of local m e across the genome, could allow for more detailed inferences about the genetic architecture of local adaptation. Our results indicate that, if local adaptation is indeed polygenic and heterogeneous across the genome, and adaptive divergence is substantial, the scope for detailed inference of the effects of individual loci contributing to local adaptation may be limited, as localized differentiation will often to a large extent reflect the genome-wide barrier effect. Nevertheless, just to what extent we can learn about the detailed genetic architecture of local adaptation and RI from observed genomic differentiation remains an open question.

The genetic architecture of local adaptation at equilibrium
When there is appreciable gene flow, only a subset of the divergently selected loci that underlie a locally adaptive trait will actually exhibit substantial differentiation at migration-selection balance, and the DFE of these loci need not be representative for the DFE associated with all loci underlying the trait (that is, the realized genetic architecture of local adaptation may differ to greater or lesser extent from the genetic architecture of locally adaptive traits (Yeaman and Whitlock 2011)). We find that when selection is fairly weak (Ls is small), the subset of divergently selected loci that exhibit significant differentiation at migrationselection balance constitutes a more biased subset than when Ls is large (fig. 5). Similarly, the DFE at migration-selection balance shifts more and more to larger selection coefficients when the strength of gene flow increases ( fig. 5), whereas the effect on the distribution of dominance coefficients depends on the correlation between s and h in the DFE underlying locally adaptative traits (figs. 5, S19, S20, S21). Importantly, the maintenance of polygenic migration-selection equilibrium itself generates correlations between selection and dominance coefficients. Correlations between selection and dominance coefficients are often discussed in the context of adaptation, with different factors influencing the relationship between the homozygous effect and dominance deviation in the DFE of new mutations (Agrawal and Whitlock 2011;Manna et al. 2011), the DFE of the standing genetic variation (Orr and Betancourt 2001;Zhang et al. 2004), and the DFE of variants fixed during adaptation (Orr 2010). We show that migration-selection balance may be another source of s-h correlation: especially when RI is low and gene flow rather strong, recessive alleles that are divergently maintained at equilibrium tend to have higher selection coefficients than dominant alleles, even when no such correlation exists a priori ( fig. 6).

Limitations of the model
Throughout, we have ignored physical linkage of the loci under selection. When recombination is strong relative to selection (i.e. linkage is weak), the above theory should work with minor modifications. However, accounting for tight linkage using an approach like ours, based on plugging in a suitable effective migration rate in single locus theory, is likely not straightforward. Associations between tightly linked loci will be broken down by recombination at rates comparable to or slower than their elimination by selection, increasing the strength of coupling (Barton 1983;Kruuk et al. 1999) and rendering the separation of time scales argument inappropriate. The relative importance of the barrier due to linked versus unlinked loci depends on the total map length as well as the number of selected loci. For organisms with a limited number of chromosomes (e.g. Drosophila), linkage may be important, whereas with larger chromosome numbers, most of the barrier may be due to unlinked loci. Secondly, we have assumed that local fitness is determined by an additive trait under directional selection, considering a history where a rapid polygenic selection response has driven allele frequencies at L loci up near fixation. Thus, we effectively assume that the locally adapted population is not yet close to a fitness optimum, so that we can ignore stabilizing selection. However, when there is abundant standing variation, a polygenic selection response may initially only involve subtle changes in allele frequencies (Sella and Barton 2019;Hayward and Sella 2022), and there may be considerable genetic redundancy (Yeaman 2015;Barghi et al. 2020), leading to a scenario that is quite different from the one assumed in this paper. The extent of reproductive isolation that can be maintained when these aspects of polygenic adaptation become important remains unclear and likely requires different approaches. More generally, our focus on the maintenance of polygenic local adaptation and the reproductive isolation it causes provides only half of the picture, as we have both ignored the initial polygenic response, and the further building up of divergence in the face of gene flow. We have considered how a given genetic architecture underlying local adaptation results in observable patterns of adaptive differentiation at equilibrium, but remain ignorant about just what sort of genetic variation is likely the source of local adaptation, neither have we considered how the resulting barrier to gene flow promotes further divergence. All these are important topics deserving further study if we are to understand how populations can remain locally adapted when subjected to maladaptive gene flow, and, ultimately, the adaptive processes that could drive the origin of new species.

Acknowledgements
This work was funded by the European Union (ERC BryoFit 101041201 granted to CF). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European S1 Supplementary figures Figure S1: Critical swamping thresholds for the multilocus model. Equilibria of the multilocus system correspond to the zeros of f (p) = hq + (1 − 2h)q 2 − m e /s. Examples for f (p) in the case with dominant local adaptation (top row), additive local adaptation (middle row) and recessive local adaptation (bottom row) near the critical point. The stable equilibrium is indicated by a filled dot, the unstable by an unfilled dot. When there is bistability, i.e. both a stable and unstable equilibrium, the critical migration rate at which the two equilibria collide and cease to exist corresponds to the value of m for which f (p) reaches its maximum in the critical point, so that both f (p) = 0 and f ′ (p) = 0 are satisfied. Figure S2: Critical equilibrium differentiation (p c , the frequency of the locally beneficial allele on the island just before swamping) and critical migration rate (m c ) for intermediate dominance (0 ≤ h ≤ 1) and low to appreciable divergence (Ls ≤ 1.5). The solid white line marks the region of parameter space where the system exhibits bistability (i.e. a sharp swamping threshold at a critical differentiation level p c > 0). The dashed lines mark h = 1/3 and h = 2/3. For h > 2/3, bistability occurs for all m. For 0 < h < 1/3, the minimum Ls for which bistable behavior is observed increases, with increasing h, after which it quickly falls.   ∈ [5, 10, 25, 50]) and h varies over rows (h ∈ [0, 0.5, 1]). We assumed k = 5 diploids per haploid individual and set N = N e /2k +N e so that the desired N e = N e s/(Ls/L) is obtained. Allele frequencies for the individual-based simulations are obtained by simulating for 110000 generations, sampling every 5 generations after discarding the first 60000, and averaging across loci. For each L we simulate n replicates so that nL = 50. The mutation rate was set to u = 0.005s.                      Where, from eq. 2, we havē Assuming the intensity of selection per generation is weak (all s are small) and that ϵ measures the generation time, we obtain a continuous-time model of allele frequency change by considering the per-generation change in allele frequency and taking the limit as ϵ goes to zero, specifically, plugging eq. 5 and eq. 4 in eq. 3, we geṫ This has the same form as the classical diploid or haploid continuous-time model of allele frequency change 2 but with marginal and mean Malthusian fitnesses given bys i ands respectively.
As an example, consider the biallelic case with alleles A 0 and A 1 , so that s 0 = s 00 = 0, p 0 = p and p 1 = q. Note that we shall always assume s ij = s ji . We have from eq. 7 s 1 = s 1 + s 01 p + s 11 q ands = s 1 q + 2s 01 pq + s 11 q 2 . Some algebra shows that we can write the ODE for the frequency of the selected allele (A 1 ) aṡ where s a = s 1 + s 01 and s b = s 11 − 2s 01 . As expected, this is the same dynamical law as for the strictly diploid model, in which case the dynamics of the selected allele are given by eq. 8 but with s a = s 01 . This enables us to identify a pair of 'effective' selection coefficients, so that, for weak selection, a diploid biallelic model with parameters s * 01 and s * 11 yields the same allele frequency dynamics 3 as a haplodiplontic model with parameters s 1 , s 01 and s 11 .

S2.2 Equilibrium structure of the mainland-island model
We describe the equilibrium structure of the haplodiplontic single-locus deterministic mainlandisland model for the biallelic case. The dynamics are given by the ODĖ = m∆q + pq(s 1 + s 01 + (s 11 − 2s 01 )q), where q is the frequency of the locally selected allele A 1 , and p = 1 − q is the frequency of the allele with relative fitness of 1 on the island when homozygous.
When m = 0 (no migration), there will be an admissible fixed point when either of the following conditions holds s 01 > −s 1 and s 01 > s 1 + s 11 (12) s 01 < −s 1 and s 01 < s 1 + s 11 , i.e. when there is ploidally antagonistic selection, diploid over-or underdominance, or both. The fixed point is obtained atp This will correspond to a stable polymorphism whenever s 01 > 0. This case was first analyzed in a discrete-time model by Scudo (1967). Now consider m > 0. We shall assume that ∆q = q mainland − q = 1 − q = p, i.e. the mainland is fixed for the locally selected allele. To describe the equilibrium behavior, it is helpful to factor the dynamical law asq Linear stability at a fixed pointq is determined by If s b = 0, we have an effectively haploid model (i.e. genic selection), and will have a stable polymorphic equilibrium atq = −m/s a whenever m < −s a , and a stable boundary equilibrium atq = 1 when m > −s a . When s b ̸ = 0, polymorphic equilibria, when they exist, will correspond to the roots of the quadratic expression in parentheses in eq. 15. These are We have the following biologically relevant equilibria: i.q = 1 (swamping) is always a stable equilibrium when m > −(s a + s b ).
ii. When 0 < m < −(s a + s b ) there is always a single stable polymorphic equilibrium at q − (dark gray zone in fig. S27) and q + will not lie in [0, 1].

Figure S27
: Equilibrium and stability behavior of the single-locus biallelic haplodiplontic mainlandisland model. The dark gray zone indicates the parameter region where there is a single protected stable polymorphic equilibrium. The light gray zone shows the parameter region where there is both a stable (but unprotected) and an unstable polymorphic equilibrium. s 1 , s 01 and s 11 are the haploid and diploid selection coefficients for the invading allele (A 1 ) on the island.
iii. When −(s a + s b ) < m < s b and 4s b /m < (s a /m) 2 ⇐⇒ 4m < s 2 a /s b , there is, besides the stable boundary equilibrium atq = 1, an unstable (repelling) equilibrium at q + , and a stable polymorphic equilibrium at q − (light gray zone in fig. S27).
The relation between the key parameters s a /m and s b /m and the equilibrium behavior of the system when m > 0 is illustrated in fig. S27.
When condition (iii) holds, sharp thresholds for swamping are observed, in which case there is a certain critical allele frequency p c below which no local adaptation cannot be maintained whatever the migration rate. We can ask for which degree of dominance such sharp thresholds for swamping can possibly be observed. From eq. 15 we see that at an equilibrium which does not correspond to p = 0, the condition f (q) = m + s a q + s b q 2 = 0 holds. A sufficient condition for observing a sharp threshold is that f obtains a maximum for some q < 1, hence that f ′ (1) > 0 where f ′ (q) = s a + 2s b q. This will be the case whenever s a + 2s b > 0. In the diploid case with the usual parameterization where s a = −sh and s b = −s(1 − 2h), this shows that critical behavior is expected as soon as h > 2/3.

S2.3 Fixed point iteration algorithm
Our approximations yield a system of equations for the expected allele frequencies and heterozygosities on the island which are coupled through the gff. To calculate expected allele frequencies and allele frequency distributions at equilibrium, we solve the system selfconsistently by performing a fixed point iteration. In words: for a given initial set of allele frequencies and heterozygosities, we calculate the gff at each locus using eq. 6; using these gff values, we next calculate expected allele frequencies and heterozygosities at each locus using numerical quadrature. This process is repeated until convergence. The algorithm is more formally outlined in algorithm S1. Algorithm S1 Fixed point iteration for calculating the expected allele frequency and expected heterozygosity on the island. e,j , s j )dp 7: (p j q j ) (n) ← 1 0 p(1 − p) ϕ(p; N e , u, m (n) e,j , s j )dp n ← n + 1 11: end while 12: return p (n) , (pq) (n)

S2.4 Swamping thresholds for the deterministic multilocus model
We now take a closer look at the equilibria of eq. 11 and their critical behavior. Clearly, p = 0 is always a solution of eq. 11, and it will correspond to a locally stable equilibrium (i.e. swamping) whenever m/s > 1 − h. Other equilibria, when they exist, are given by the zeros of the function f (p) = hq + (1 − 2h)q 2 − m s g [p] (17) Note that g[p] > 0 and we assume s > 0, so that for any fixed h, as m increases, there will indeed be a critical migration rate beyond which f (p) < 0, from which point onwards the only stable equilibrium will be p = 0. At the critical point, the equilibrium allele frequency will satisfy the additional constraint f ′ (p) = 0 (see fig. S1), i.e.
f ′ (p) = h + 2(1 − 2h)q − 2Lm(1 − 3h − 2(1 − 2h)q)g[p] = 0 We can solve eq. 17 for g [p], and then plug in g [p] in eq. 18. This yields a cubic polynomial in p which can be solved for the allele frequency p c at the critical point: We can then plug p c into eq. 17 and solve for m c /s. While the general expressions yield not much insight, we can focus on a number of special cases.
Firstly, in the additive case (h = 0.5), a critical point different from 1 − h appears when Ls > 1, in which case the equilibrium frequency at the critical point will be p c = 1 − 1/Ls. The corresponding critical migration rate is m c s = e Ls−1 2Ls . Figure S28: Different apparent equilibria depending on initial conditions. In the left plot, the yellow line (p 0 = 1) indicates the expected allele frequencies as determined using the fixed point iteration of algorithm S1, starting with p (0) = (1, 1, . . . , 1) (i.e. secondary contact, maximal initial differentiation), whereas the green line assumes p (0) = (0, 0, . . . , 0) (no initial differentiation). The dots show results from individual-based simulations with the same initial conditions (50000 generation, keeping the last 25000 and subsampling every 5 generations). The three plots on the right show the evolution of the fixed point iteration for different initial initial conditions p (0) = (p 0 , p 0 , . . . , p 0 ) for three values of m. Note the bifurcation of the dynamical system defined by the algorithm: for m/s = 0.2 and m/s = 0.4 there is a single globally stable fixed point, whereas for m/s = 0.3, there are two locally stable fixed points.

S2.6.2 Logistic model
In the logistic model, we assume, for i = 1, . . . , L, logit h i |s i ∼ Normal(a + b log s i , σ 2 ), where logit h = log h 1−h is the logit transform. In other words, we assume h to be distributed according to a linear regression on log s with slope a and intercept b, on a logit scale. The where N(·; µ, σ) and G(·; κ, λ) denote the density functions for the Normal distribution (with mean µ and standard deviation σ) and Gamma distribution respectively. Instead of setting the a and b parameters directly, we parameterize the regression by choosing two reference pointst, s (1) and s (2) , together with their respective expected dominance coefficients h (1) = E[h|s (1) ] and h (2) = E[h|s (2) ] using

S2.6.3 Model after Caballero and Keightley (1994)
In the model of Caballero and Keightley (1994) (see also Zhang et al. (2004) and discussion in Agrawal and Whitlock (2011)), referred to as CK94, we assume It should be noted that this distribution is supposed to be a reasonable model for dominance coefficients of deleterious mutations at mutation-stabilizing selection equilibrium, where mutations of large effect segregating at appreciable frequencies tend to be recessive. In our case, we regard this as the distribution of dominance coefficients of mutant alleles on the mainland that constitutes the standing variation which is the source of locally adaptive alleles during the initial polygenic response (which we do not explicitly model). These are hence the dominance coefficients of the locally beneficial alleles on the island (assuming dominance coefficients to be constant across environments and genetic backgrounds). The h i as we defined them are however the dominance coefficients of the invading wild-type alleles from the mainland over the locally beneficial ones, so that we use h i = 1 − h * i where the h * i are distributed according to the CK94 model. We set the K parameter so that E[h * ] =h for someh, i.e. K = λ (2h) − 1 κ − 1 .
The marginal density for h * is where γ is the lower incomplete gamma function. This model is also illustrated in fig. S29 (C) and fig. S30. To study the effect of a negative correlation between s and h (i.e. where strongly selected locally beneficial alleles tend to be dominant), we use the same model, but with h i = h * i . We refer to this model as CK94 * (see fig. S26).