Establishment of local adaptation in partly self-fertilizing populations

Populations often inhabit multiple ecological patches and thus experience divergent selection, which can lead to local adaptation if migration is not strong enough to swamp locally adapted alleles. Conditions for the establishment of a locally advantageous allele have been studied in randomly mating populations. However, many species reproduce, at least partially, through self-fertilization, and how selfing affects local adaptation remains unclear and debated. Using a two-patch branching process formalism, we obtained a closed-form approximation under weak selection for the probability of establishment of a locally advantageous allele (P) for arbitrary selfing rate and dominance level, where selection is allowed to act on viability or fecundity, and migration can occur via seed or pollen dispersal. This solution is compared to diffusion approximation and used to investigate the consequences of a shift in a mating system on P, and the establishment of protected polymorphism. We find that selfing can either increase or decrease P, depending on the patterns of dominance in the two patches, and has conflicting effects on local adaptation. Globally, selfing favors local adaptation when locally advantageous alleles are (partially) recessive, when selection between patches is asymmetrical and when migration occurs through pollen rather than seed dispersal. These results establish a rigorous theoretical background to study heterogeneous selection and local adaptation in partially selfing species.


Deriving the difference equation for allele frequency
Using the preceding definitions, one can derive the difference equations for allele frequencies 139 in each patch. Let X i and Y i be the frequencies of genotypes AA and Aa, and x i be the 140 frequency of allele A in i th patch. Our system has two degrees of freedom, given that there 141 are three genotype frequencies and they must add up to unity. Thus, we keep track of 142 frequencies of genotypes AA and Aa, while the frequency of aa genotype is by definition 143 1 − X i − Y i . Adult genotype frequencies in generation t in i th patch are: Adult genotypes may produce different number of male and female gametophytes. Let W ♂ i,k 145 and W ♀ i,k be the number of male and female gametophytes that an adult with genotype k 146 (AA, Aa and aa) produces in patch i. Thus, the frequency of female and male gametophyte 147 A in i th patch are, respectively: After pollen migration, the frequency of allele A in pollen in i th patch is: where m ij is the fraction of pollen in i th patch that come from j th patch. Pollen dispersal is 151 followed by mating, which yields offspring genotypes: Seed migration changes the genotype frequency to: where M ij is the fraction of seed in patch i originating from patch j. Genotype j in patch 154 i has probability V i,j to survive to maturity, so the frequency of allele A in the generation 155 t + 1 is: dynamics for an arbitrary mode of selection and migration, it is difficult to analyze, so we 180 focus on special population genetic cases outlined above. Each of these scenarios will have 181 an associated Jacobian J {k} : 182 x 1 (t + 1) Terms ∆V 0,i , ∆W These expressions can be made more intuitive by parameterizing fitnesses relative to geno- are the reproductive advantage of AA homozygote relative to 195 aa homozygote when selection acts on viability, female, and male fecundity, respectively.
These forms capture the fact that selfing increases the effective dominance of invading through seed production (female fitness). Taken together: The fitness of the mutant always take the same form as in a haploid model with an effective 211 selective advantages =hs.

213
The preceding approach can be readily extended to the selective scenarios where migration 214 occurs solely via pollen. One can still use the Jacobian in (7), albeit with the modified 215 migration rate to account for the facts that (a) pollen is haploid and thus do not contribute 216 in the same way as diploid seed to the gene pool of the next generation, and (b) that in all 217 three cases migration operates prior to selection rather than after: Selection terms are still parameterized according to equations (10a-10c), but effective mi-219 gration is given by: Viability effective migration rate (12a) can be intuitively understood by noting that 1 − S 221 corresponds to the fraction of pollen that is exchanged between the patches, m ij /2 corrects 222 for the fact that haploid phase contributes two-fold less alleles to the gene pool of the 223 next generation relative to the diploid progeny dispersal, and the relative fitness ratio (1 + Overall, a simple parameterization procedure allows the exploration of a wide breadth of 232 population genetic scenarios. By analyzing the stability of system 7, one can only examine conditions necessary for the 235 invading allele to escape extinction. However, we also want to know the probability of the 236 escape. We pursue this by representing the change in the number of invading allele as 237 a multi-type branching process. Selection is thus assumed to be stronger than drift (i.e., As we can assume that the number of offspring that types i and j leave are uncorrelated, we have: The first term in equation (15) recalling that M i,i = 1 − M i,j with j = i. Note that fecundity selection occurs in patch i, 280 before seed migration, whereas viability selection occurs in patch j, after seed migration.

281
However, the different orderings of migration and selection yield the same results (see Math-282 ematica notebook). We also retrieve that the mean of the four distributions are given by the 283 Jacobian terms obtained above, M i,j ∆W i,j and the variance is also inflated by 1 + F . For 284 pollen migration the form is slightly different because migrants can only contribute to the 285 next generation through outcrossing and male gametes (whereas selfed seeds can migrate).

286
The form of the PGF is thus different for the resident and migrant contribution. We have: which increase the proportion of resident offspring contributed by selfing, hence the variance.

294
Note, however, that the difference in PGFs between the two migration modes only affects 295 high-order terms of selection and migration so has no effect on branching process approxi-296 mation results, which can all be expressed with the same form using appropriate effective 297 parameters as defined in previous section (see Mathematica notebook).

299
Previous section also showed that ∆W 0,i always has the form of 1+s i , wheres i is the effective where µ i and σ i are drift and diffusion coefficients for patch i, and are given by: Note thats corresponds to effective selective advantage of mutant, andM ij is the effective 319 migration rate; These were derived in section 2.2. Diffusion coefficient is inflated by the 320 factor of 1 + F as described previously. Unfortunately, we could not find neither the explicit 321 nor approximate solution to this partial differential equation. Following treatment in Barton

322
(1987), the problem is simplified by focusing on the first phase of fixation when mutant allele 323 is rare. Thus, one can approximate coefficients as: which leads to the simplified KBE: Because the invading allele is rare, the mutants in patch i go extinct independently of one 326 another with probability e −l i , which means that the establishment probability takes the form: Substituting eqn. 22 into 21 results in an algebraic equation which captures the establishment 328 probability for an arbitrary number of initial mutants across two patches. We are interested 329 in two special cases, when mutant either starts in favored (x 1 = 1/(2N 1 ) and x 2 = 0) or 330 disfavored patch (x 1 = 0 and x 2 = 1/(2N 2 )). With this parameterization, we obtain the 331 system of two non-linear equations: By solving these for l 1 and l 2 , taking the real root, and substituting into eqn. 22, one retrieves 1, as is the case with branching approximation. As will be shown later 336 in the text, the comparison of two solutions reveals that diffusion approximation performs 337 identically when the allele starts in the favored patch and slightly better when it originates 338 in the disfavored one. Because branching process approximation is easier to work with (e.g., 339 series expansion readily yields intuitive special cases), we use this solution throughout our 340 analysis and note that the diffusion approach yields marginally better precision.
The denominator is the relative reproductive advantage of invading mutant against the of allele A in either patch 1 or 2 is given by: where ∆ =h 1 s 1 −h 2 s 2 and σ = √ ∆ 2 + M 2 withh as defined above for viability, female  that we also accounted for the reduction in P due to increased variance of offspring number.

376
Conversely, in the limit of full migration we obtain P 1 = P 2 =h 1 s 1 +h 2 s 2

1+F
, which is equivalent 377 to Nagylaki (1980), with the appropriate rescaling for dominance and partial selfing. The . Thus, unlike in a single population a locally advantageous 384 allele that is codominant in both patches is more likely to escape extinction in selfers than 385 in outcrossers. The relation between F and P 2 is more complicated, and the allele can either 386 be more likely, for high migration rates, or less likely, for low migration rates, to escape

Comparison with simulations
In general, our approximation for viability shows good correspondence with simulations 391 across different selfing rates ( Figure 1). This also holds true across different dominance co- P 1 is a monotonically decreasing function in M , because strong migration causes the spread-403 ing of the locally advantageous allele to become swamped by the resident from the disfa-404 vored patch. P 2 is non-monotonic function in M . Increasing migration initially increases the 405 probability that an invading allele escapes from the disfavored patch before it goes extinct.

406
However, at large migration rates, the spreading allele escapes the disfavored patch but is 407 swamped due to a large influx of deleterious residents from the disfavored patch. For male 408 fecundity selection, the analytical solutions for P 1 and P 2 are slightly worse against the sim-409 ulations ( Figure A11) than those provided for viability selection. Note, however that as the 410 selfing rate tends towards 1, the effective selection coefficient on male function vanishes so 411 that the strong selection condition for branching process approximation is not met anymore.  We also varied F continuously to obtain a fine-grained view of the analytic performance.
Consider decomposing the effects of selfing on the establishment of an allele affecting viability.

438
For example, the indicator β 1 (h • 1 ) tells us whether a shift to selfing increases P 1 , if selfing . Lines in the right panel are as follows. Solid black line: . Parameters for the left panel: For   The above analysis showed how selfing affects the critical migration rate, hence the conditions 570 for local adaptation. However, when local adaptation is possible, selfing also affects the 571 probability of establishment and maintenance of polymorphism. This can be analyzed by 572 quantifying how selfing affects P (A) and P (a) , which is given by: Four outcomes are possible upon the introduction of selfing. First, P (A) increases, and P (a) 574 decreases, meaning that allele A is more likely to fix across both patches in selfers than 575 in outcrossers. Second and conversely, P (A) decreases, but P (a) increases, implying that 576 selfers are less likely to fix A compared to an outcrossing population. Third, both P (A) 577 and P (a) increase upon shift to selfing, meaning that selfing increases the probability that One can intuitively understand these outcomes by noting that dominance reverses when the 583 alternative allele invades. So, if invader A is dominant, then invader a will be recessive.

584
The dominant allele A is likely to escape extinction when invading but unlikely to result 585 in protected polymorphism because allele a is recessive and thus likely to go extinct. This We assume that the number of exported pollen and the different numbers of seeds follow 783 independent Poisson distributions. The different means depends on the fecundity of the 784 parent, the viability of the seed produced and the proportion of each category: The factor 1/2 corresponds to the fact that the focal A is chosen with probability 1/2. Put 786 another way, on average each individual transmits two alleles so the contribution for a single 787 allele must be halved. Finally, we use the property that the PGF of 2X is f (z 2 ) where f (z) 788 is the PGF of X and that the PGF of the sum of independent variables is the product of the 789 PGFs. Noting g(λ, z) = e −(1−z)λ , the PGF of a Poisson distribution with mean λ, we have: which yields equation (15) of the main text.

791
From the properties of PGFs we can easily obtain the moments of the distribution for the 792 different forms of selection (see Mathematica notebook): So we retrieved the well known result that selfing increases variance in allele offspring number 794 by 1 + F . Interestingly, the full distribution presents a peculiar non-monotonic behavior for 795 high selfing, with an excess (resp. a deficit) in even (resp. odd) numbers. Note also that the Note that in subsequent analyses the order of migration and selection terms yields the same whereas for offspring contributing to the other patch: and for migrant offspring: (1 +h j s j ) for viability fecundity selection (A10a) = 0 for female fecundity selection (A10b) whereM i,j ≈ m i,j (1−S)/2 are effective migration rates as defined in the main text (equations 819 12a to 12c) and ω is a correcting factor defined in equation (13). Compared to seed migration, 820 the variance is not uniformly increased by 1 + F . As migrant offspring can only be produced 821 through outcrossing, the distribution is simply Poisson and the variance equal to the mean.

822
On the contrary, because the proportion of outcrossed offspring contributing to the resident 823 patch is reduced due to pollen migration, the variance is inflated by more than 1 + F . Eigenvectors are normalized such that: We choose a parameter in the model, such that the process is slightly supercritical when As → 0, ρ approaches unity and the second term in B (equation (A12)) can be neglected.
Recall thats i is the advantage of an allele in patch i accounting for the effect of self-840 fertilization. Assuming weak selective advantage,s 1 is taken as . All other parameters are 841 expressed in terms ofs 1 : The leading eigenvalue of M can be written as ρ = 1+cs 1 +o(s 2 1 ). Dropping the higher-order 843 terms ins 1 , we retrieve the expression forc:  If the selective disadvantage in the disfavored patch is too large or migration is too strong, 866 the spreading locally advantageous allele can be swamped by its deleterious counterpart.

867
The range of parameters that are necessary but not sufficient for a successful invasion are for this equilibrium to be locally unstable. If migration occurs prior to selection -which is 870 the case when seed disperses and selection acts on viability, then Jacobian J is: Symbol W is a place-holder for any of the three selection modes (V , W ♀ , or W ♂ ). If 872 selection occurs prior to migration -which happens when seed disperses and selection acts 873 on sexual components-then Jacobian is: However, the eigenvalues of these two matrices are identical, so we here use (A18). The 875 equilibrium is locally unstable whenever the leading eigenvalue of J is greater than unity: (because A is advantageous in the first patch), and 0 < ∆W 0,2 < 1 (because A is deleterious 878 in the second patch). It is important to note that these inequalities hold regardless of the 879 mode of selection. More formally: Conversely to the previous case, 0 < ∆W 1,1 < 1 (as invading allele a is deleterious in the 887 first patch) and ∆W 1,2 > 1 (given that a is beneficial in the second patch). Once the allele A has escaped extinction, it can either fix in both patches or be maintained for finite number the same inequalities hold as in the case of seed dispersal. However, under male fecundity 902 selection migration rates have to be: for both m 12 and m 21 , or:      and disfavored dominance (β k (h 1 ) and β k (h 2 )), and the effective migration rate (β k (M )). • β k (h 2 ): Parameterize migration according to eqns. 12a-12c, and then set F = 0 to Figure A17: The consequences of a shift to selfing on establishment of local adaptation under female fecundity selection and seed dispersal. Color-coding and parameters as in the main text and parameters as in the main text. Figure A18: The consequences of a shift to selfing on establishment of local adaptation under male fecundity selection and seed dispersal. Color-coding and parameters as in the main text and parameters as in the main text.  For the sake of clarity, Table 1  The effective selective advantage of invading allele A with selection via k th fitness component