Species often encounter, and adapt to, many patches of similar environmental conditions across their range. Such adaptation can occur through convergent evolution if different alleles arise in different patches, or through the spread of shared alleles by migration acting to synchronize adaptation across the species. The tension between the two reflects the constraint imposed on evolution by the underlying genetic architecture versus how effectively selection and geographic isolation act to inhibit the geographic spread of locally adapted alleles. This paper studies the balance between these two routes to adaptation in a model of continuous environments with patchy selection pressures. We address the following questions: How long does it take for a novel allele to appear in a patch where it is locally adapted through mutation? Or, through migration from another, already adapted patch? Which is more likely to occur, as a function of distance between the patches? What population genetic signal is left by the spread of migrant alleles? To answer these questions we examine the family structure underlying migration–selection equilibrium surrounding an already adapted patch, treating those rare families that reach new patches as spatial branching processes. A main result is that patches further apart than a critical distance will likely evolve independent locally adapted alleles; this distance is proportional to the spatial scale of selection (σ / s m, where σ is the dispersal distance and sm is the selective disadvantage of these alleles between patches), and depends linearly on log (sm/μ), where μ is the mutation rate. This provides a way to understand the role of geographic separation between patches in promoting convergent adaptation and the genomic signals it leaves behind. We illustrate these ideas using the convergent evolution of cryptic coloration in the rock pocket mouse, Chaetodipus intermedius, as an empirical example.
The convergent evolution of similar phenotypes in response to similar selection pressures is a testament to the power that selection has to sculpt phenotypic variation. In some cases this convergence extends to the molecular level, with disparate taxa converging to the same phenotype through parallel genetic changes in the same pathway, genes, or even by precisely the same genetic changes [1–3]. Convergent adaptation also occurs within species, if different individuals adapt to the same environment through different genetic changes. There are a growing number of examples of this in a range of well studied organisms and phenotypes [4]. Such evolution of convergent phenotypes is favored by higher mutational input, i. e. , higher total mutational rate and/or population size [5]. The geographic distribution of populations can also affect the probability of parallel mutation within a species: a widespread species is more likely to adapt by multiple, parallel mutations if dispersal is geographically limited, since subpopulations will adapt via new mutations before the adaptive allele arrives via migration [6]. Standing variation for a trait can also increase the probability of convergence, as this increases the probability that the selective sweep will be soft (beginning from a base of multiple copies), which leads to genetic patterns similar to convergent adaptation [7,8] whether or not the copies derive from the same mutation, Intuitively, convergence is also more likely when geographically separated populations adapt to ecologically similar conditions. The probability that convergent adaptations arise independently before adaptations spread between the populations by migration will be larger if these adaptive alleles are maladapted in intervening environments, since such adverse conditions can strongly restrict the spread of locally adapted alleles [9]. One elegant set of such examples is provided by the assortment of plant species that have repeatedly adapted to patches of soil with high concentrations of heavy metals (e. g. , serpentine outcrops and mine tailings) [10–13]; the alleles conferring heavy metal tolerance are thought to be locally restricted because they incur a cost off of these patches. Similarly, across the American Southwest, a variety of species of animals have developed locally adaptive cryptic coloration to particular substrates, e. g. , dark rock outcrops or white sand dunes [14]. One of the best-known examples is the rock pocket mouse (Chaetodipus intermedius), which on a number of black lava flows has much darker pelage than on intervening areas of lighter rock [15]. Strong predator-mediated selection appears to favour such crypsis [16], and, perhaps as a result of this strong selection against migrants, at least two distinct genetic changes are responsible for the dark pigmentation found on different outcrops [17]. Similar situations have been demonstrated in other small rodent systems [18–20] and in lizards [21]. In this paper, we study this situation, namely, when a set of alleles provide an adaptive benefit in geographically localized patches, separated by inhabited regions where the alleles are deleterious. The main questions are: Under what conditions is it likely that each patch adapts in parallel, i. e. , convergently through new mutation, and when is it likely that migration carries these alleles between patches? How easy will it be to detect adaptive alleles that are shared by migration, i. e. , over what genomic scale will a population genetic signal be visible? We work in a model of continuous geography, using a variety of known results and new methods. In the section Establishment of a locally adaptive allele due to mutational influx we develop a simple approximation, Eq (2), for the rate at which patches become established by new mutations. The most important conceptual work of the paper occurs in the section Establishment of a locally adaptive allele due to migrational influx, where we develop an understanding of the process by which alleles move from an existing patch to colonize a novel patch, culminating in Eq (11) for the rate at which this happens. We combine these two results in the section The probability of parallel adaptation between patches to discuss the probability of parallel adaptation, Eq (12). To understand the genomic signal left by adaptations shared by migration, in the section Length of the hitchhiking haplotype, we study the time it will take for an allele to transit between patches, (Eq (18) ), and thus the length of haplotype that hitchhikes with it (Eq (19) ). Finally, in the section Applications we apply this work to understand the geographic scale of convergent evolution in Chaetodipus intermedius.
Consider a population spread across a landscape to which it is generally well adapted, but within which are patches of habitat to which individuals are (initially) poorly adapted. (When we refer to “patches” it is to these pieces of poor habitat.) Suppose it takes only a single mutational change to create an allele (B) that adapts an individual to the poor habitat type. The required change could be as specific as a single base change, or as general as a knockout of a gene or one of a set of genes. This change occurs at a (low) rate of μ per chromosome per generation, and has fitness advantage sp relative to the unmutated type (b) in these “poor” habitat patches. Outside of these patches the new allele has a fitness disadvantage of sm in the intervening areas, with sp and sm both positive. (Here we define “fitness” to be the intrinsic growth rate of the type when rare.) We assume that the disadvantage sm is sufficiently large that, on the relevant timescale, the allele is very unlikely to fix locally in these intervening regions. (The case where sm = 0 requires a different approach, which we do not treat here.) We are interested in the establishment of mutations in the “poor” patches by either migration or mutation, which depends on whether the allele can escape initial loss by drift when rare. Therefore, we need not specify the fitness of the homozygote; only that the dynamics of the allele when rare are determined by the fitness of the heterozygote. More general dominance will only make small corrections to the dynamics until initial fixation, with the exception of the recessive case, which we omit. In other words, we follow the literature in treating the diploid model as essentially haploid. We also assume population density ρ is constant across the range (even in the “poor” patches). The variance in offspring number is ξ2, and that the mean squared distance between parent and child is σ2 (i. e. , σ is the dispersal distance). We will deal with migration by immediately appealing to the central limit theorem, treating the total distance traveled across t dispersal events as Gaussian with variance tσ2, and do not discuss the (interesting) cases where rare, long-distance dispersal events are more important (see e. g. , [6,22,23] for discussion). We first compute the time scale on which a new B mutations will appear and fix in a single, isolated poor habitat patch in which no B allele has yet become established. As we are interested in patches where local adaptation can occur, we will assume that the patch is larger than the cutoff for local establishment mentioned above. Let p (x) be the probability that a new mutant B allele that arises at location x relative to the center of the patch fixes within the poor habitat patch. Under various assumptions, precise expressions can be found for p (x) [29], but results will be more interpretable if we proceed with a simple approximation. The total influx per generation of mutations that successfully establish is the product of population density ρ, mutation rate μ, and the integral of p (x) over the entire species range: λ mut = ρ μ ∫ p (x) d x. (1) This depends in a complicated way on the patch geometry and selection coefficients, but still scales linearly with the mutational influx density ρμ. If we consider a patch of area A, whose width is large relative to σ / 2 s m, then a reasonable approximation is to ignore migration, so that p (x) = pe ≈ 2sp/ξ2 within the patch, and p (x) = 0 otherwise. This approximates the integrand p (x) in Eq (1) by a step function, which will be a good approximation if the patch is large relative the scale over which p (x) goes from 0 to pe, or if p (x) is close to pe at some point and is symmetric about the edge of the patch. We examine this more generally via exact calculation of p (x) in the section Numerical calculation of the probability of establishment. The rate at which mutations arise and colonize a patch of area A is then λ mut ≈ 2 s p ρ A μ / ξ 2, (2) i. e. , the product of mutational influx in the patch (ρAμ) and the probability of establishment (pe ≈ 2sp/ξ2). If this rate is low, then the time (in generations) until a mutation arises that will become locally established within the patch is exponentially distributed with mean 1/λmut. Assuming that once a mutation becomes established it quickly reaches its equilibrium frequency across the patch, the time scale on which new patches become colonized by the B allele from new mutation is therefore approximately 1/λmut. Now suppose that there are two patches separated by distance R (i. e. , the shortest distance between the two is R). If the B allele has arisen and become established in the first patch, but has not yet appeared in the second, we would like to know the time scale on which copies of B move between patches by migration (supposing that no B allele arises independently by mutation in the meantime). To determine this, we study the fine-scale genealogy of alleles that transit between patches, obtaining along the way other useful information about the genealogy of B alleles. Doing this we arrive at Eq (11) for the rate at which an allele established in one patch spreads to a neighboring one by migration. Migration–selection balance ensures that there will be some B alleles present in the regions where they are disadvantageous, but only rarely, far away from the patch where B is established. Denote the expected frequency of allele B at displacement x relative to the patch by q (x), and assume that the shape of the patch is at least roughly circular. Following [24] or [9], one can write down a differential equation to which q (x) is the solution, and show that q (x) decays exponentially for large ∣x∣, with a square-root correction in two dimensions: q (x) ≈ C (| x | 2 s m / σ) − (d − 1) / 2 exp (− | x | 2 s m / σ) for large | x |, (3) where d is the dimension (d = 1 or 2), and C is a constant depending on the geographic shape of the populations and the selection coefficients. In applications we fit C to the data; to get concrete numbers from theory we take C = 1 if necessary. As J. Hermisson pointed out in comments on an earlier draft, this functional form has a simple intuition: local migration leads to the exponential decay, since if each migrant at distance x produces an average of c descendants that make it to distance x + 1, then in one dimension, the number of migrants must decay as exp (−cx); and the square-root term in two dimensions enters because the available area grows with x. The “renewal” aspect of this same argument suggests that for Eq (3) to hold, ∣x∣ must be larger than a few multiples of σ. These asymptotics fit simulations quite well, as shown in Fig 1. To be clear about the assumptions implicit here, below we provide a derivation of Eq (3) in section The equilibrium frequency, as well as justification for the asymptotics below in section Asymptotic solution for the equilibrium frequency. In one dimension, the equation can be solved to give q (x) in terms of Jacobi elliptic functions [37]; see the Supporting Information. This expected frequency q (x) is the time-averaged occupation frequency, i. e. , the total number of B alleles found across T generations near location x, per unit area, divided by T. If, for instance, q (x) =. 01 and the population density is ρ = 10 individuals per unit area, then in a survey tract of unit area around x we expect to find one individual every 10 generations—or, perhaps more likely, 10 individuals every 100 generations. This points to an important fact: close to the original patch, the “equilibrium frequency” q (x) describes well the density of the B allele at most times, but far away from the patch, the equilibrium is composed of rare excursions of families of B alleles, interspersed by long periods of absence. An example of one of these rare excursions is shown in Fig 2. The relevant time scale on which B alleles migrate between patches is given by the rate of appearance of such families. This suggests decomposing the genealogy of B alleles into families in the following way. First, consider the genealogy of all B alleles that were alive at any time outside of the original patch. This is a collection of trees, each rooted at an allele living outside the patch whose parent was born inside the patch. Next, fix some intermediate distance r0 from the established patch, and erase from the genealogy every allele that has no ancestors living further away than r0 to the patch. This fragments the original set of trees into a collection of smaller trees that relate to each other all B alleles living outside of r0 at any time, and some living inside of r0; we refer to these trees as “families”. If r0 is large enough that there are not too many families and interactions between family members are weak, then these families can be approximated by a collection of independent spatial branching processes whose members are ignored if they enter the original patch, illustrated in Fig 3. (This can be made formal in the limit of large population density, also taking r0 large enough that the chance of reaching the original patch is small.) The opportunity for adaptation depends on how often these families of B alleles encounter the new patch. Suppose that the area occupied by the new patch is S. We can learn about the rate of arrival of families at S from the time-averaged number of B alleles expected in S were it not a patch (i. e. , if B was still maladaptive in S), which from Eq (3) is ρ q (S) = (outflux of families) × (mean occupation in new patch per family), (4) where q (S) = ∫S q (y) dy. The quantity we wish to compute, the effective rate at which families of migrant B alleles establish in the new patch, is λ mig (S) = (outflux of families) × (probability a family establishes in new patch). (5) The quantities on the right-hand side depend implicitly on r0, but their product does not. Now that we have expressions for the mean rates of adaptation by new mutation, Eq (2), and by migration from an already colonized patch, Eq (11), it seems helpful to step back and review the assumptions underlying the asymptotic results we have used, or will use below. Our results should hold exactly in the limit of large, circular patches sufficiently far apart, large population density, and small selection coefficients of equal magnitude. To be more precise, the mutation rate results most obviously apply if sm ≈ sp, if the patch diameter is large relative to σ / min (s m, s p), and the patch is not too far from circular. If sm ≪ sp then a strip of width σ / s m should be added to the outside of the patch in computing A for Eq (2), and if sm ≫ sp, a strip of width σ / s p should be subtracted from the patch. As for the migration rate, we assume that each patch is large enough to support a stable population of B alleles (A 1 / d > (σ / 2 s p) tan - 1 (s m / s p) ). The geometry and size of the patch will also affect the approximation of Eq (48). Next, in using Eq (3), we assume that the inter-patch distance is large relative to the characteristic length (R > σ / s m), and that local drift is not so strong that the equilibrium frequency is at least approximately attained (using Wright’s effective neighborhood size, sm ≫ 1/ (ρσd) ). The last requirement is necessary because if the B allele fixes in large neighborhoods where it is deleterious, we cannot approximate its dynamics via a branching process. We also neglect the time for migration-selection equilibrium to be reached. As discussed above, we also assume that migration to a new patch takes place over a number of generations; if there are sufficient rare, long-distance migration events that would move between patches in a single hop, this would require a different analysis. To test the robustness of our results to the various approximations we used, we used individual-based simulations on a one-dimensional lattice of 501 demes, with ρ haploid individuals per deme, dispersal to nearby demes with σ = 0. 95, and run for 25,000 generations. More details are given below in section Simulation methods, and the number of simulations used and parameter values are given in supplementary S1 and S2 Tables. To estimate the mean time until adaptation by mutation, we used one centrally located patch of 99 demes and a mutation rate of μ = 10−5, and to estimate the mean time until adaptation by migration, we used two centrally located patches of 99 demes separated by a varying number of demes, with one patch initialized to have the B allele at frequency 0. 8. In each case, we measured the “time to adaptation” as the amount of time until at least 100 B alleles were present in the previously unadapted patch. Fig 4 summarizes how the results compare to theory, excluding parameter combinations that violate the assumptions discussed above, or where a majority of simulations did not adapt by 25,000 generations. S1 and S2 Figs show all times, and depictions of typical simulation runs are shown in S3, S4, S5, S6, S7, S8, S9 and S10 Figs. The agreement is reasonable for both. The theoretical value 1/λmut underestimates the mean time to mutation by a factor of around 2 that increases slightly with sm. This is to be expected for two reasons: First, we compute the time to reach 100 B alleles, while theory predicts the time until the progenitor of those 100 B alleles arose. Second, the expression for λmut neglects effects near the boundary of the patch, and the larger sm is, the harder it is for mutations that arise near the edge of the patch to establish. The theoretical values 1/λmig again has a margin of error of a factor of about 2. S10 Fig shows simulations at more parameter values. We now turn to the question of whether the separated patches adapt by parallel genetic changes or by the spread of a migrant allele between the patches. As only a single mutation is required for individuals to adapt to the poor habitat patch, subsequent mutations that arise after an allele becomes established in the patch gain no selective benefit. Similarly, an allele introduced into a patch by migration will not be favored by selection to spread, if the patch has already been colonized. Therefore, mutations selectively exclude each other from patches, over short time scales, and will only slowly mix together over longer time scales by drift and migration, an assumption we also made in [6]. In reality, different mutations will likely not be truly selectively equivalent, in which case this mixing occurs over a time-scale dictated by the difference in selection coefficients. We assume that once a B allele is introduced (by migration or mutation) it becomes established in the poor habitat patch rapidly if it escapes loss by demographic stochasticity. Under this approximation, and the “selective exclusion” just discussed, after one of the patches becomes colonized by a single B allele, the other patch will become colonized by migration or mutation, but not by both. As such, the question of how the second patch adapts simply comes down to whether a new mutation or a migrant allele is the first to become established in the second patch. To work out how adaptation is likely to proceed, we can compare Expressions (2) and (11) above for the effective rates of adaptation by new mutation and by migration. We work in one dimension, as the square root term appearing for two dimensions is relatively unimportant. We first consider the order of magnitude that our parameters need to be on in order for adaptation via mutation or migration to dominate. Let γ = min (1, sm/pe) and w = A/A′—since A is the area of the not-yet-adapted patch and A′ is the area of its closest σ / 2 s m strip, w is approximately the width of the patch in units of σ / 2 s m. Effective migration and mutation rates will be on the same order if A μ ≈ 2 A ′ γ s m exp (- R 2 s m / σ), where R is the distance between the patches. In other words, the migrational analogue of “mutational influx” (μρA) is 2 ρ A ′ γ s m exp (- R 2 s m / σ), which depends fairly strongly on the selective disadvantage sm between patches. Equivalently, the rates are roughly equal if R / σ = log (2 γ s m / (w μ) ) / 2 s m, which gives the critical gap size past which adaptation will be mostly parallel in terms of selection coefficients, patch width, and mutation rate. If we take μ = 10−5, the patch width to be ten times the length scale σ / 2 s m so w = 10 (and −log (wμ) ≈ 9. 2), and sm > pe so that γ = 1, then adaptation is mostly parallel for patches separated by gaps larger than R / σ > (9. 2 + log (2 s m) ) / 2 s m. If the selective pressure is fairly strong (say, sm =. 05), then for convergence to be likely, the distance between patches must be at least 21. 8 times the dispersal distance (i. e. , R/σ > 21. 8). If sm is much smaller, say sm =. 001, the required distance increases to R/σ > 67. If the mutation rate is higher—say, μ = 10−3—the required separation between patches is only reduced to R/σ > 7. 3 with sm =. 05. If sm =. 001, then with this higher mutational influx this calculation predicts that mutation will always be faster than migration—however, this should be taken with caution, since as discussed above, this model does not hold if R is of the same order as σ or if sm is small enough that local drift is more important. We can go beyond these rough calculations to find the probability of parallel adaptation if we are willing to take our approximations at face value. Doing so, the time until the first of the two patches is colonized by the B allele will be approximately exponentially distributed with mean 1/ (2λmut). Following this, the time until the second patch is subsequently colonized (via either migration or new mutation) will be exponentially distributed with mean 1/ (λmut + λmig). Both these rates scale linearly with population density (ρ) and the probability of establishment of a rare allele (pe ≈ 2sb/ξ2, for pe > sm), so that the overall time to adaptation is robustly predicted to increase with offspring variance (ξ2) and decrease with population density and beneficial selection coefficient. Furthermore, the probability that the second adaptation is a new mutation, i. e. , the probability of parallel adaptation, with γ = min (1, sm/pe), is λ mut λ mut + λ mig (R) = w μ / (2 s m) w μ / (2 s m) + C γ (R 2 s m / σ) − (d − 1) / 2 exp (− 2 s m R / σ), (12) so that the probability of parallel mutation should increase approximately logistically with the distance between the patches, on the spatial scale σ / 2 s m. We tested this using the same simulations as Fig 4, by using the empirical distributions of the respective times to adaptation to estimate the probability, for each parameter combination, that a new, successful mutation appears before a successful migrant arrives from another, already adapted patch. The results are compared to Eq (12) in Fig 5. If a patch adapts through new mutation or a rare migrant lineage, the genomic region surrounding the selected site will hitchhike along with it [42], so that at least initially, all adapted individuals within a patch share a fairly long stretch of haplotype. Pairs of adapted individuals within a patch will initially share a haplotype of mean genetic length of about log (ρAsp) /sp around the selected site, as long as the patch is reasonably well mixed by dispersal (otherwise see [43]). This association gets slowly whittled down by recombination during interbreeding at the edge of the patch, but there will always be longer LD nearby to the selected site [44]. When an already adapted patch colonizes another through migration, the newly colonized patch will inherit a long piece of haplotype around the selected site from the originating patch. The genetic length of this haplotype is roughly inversely proportional to the number of generations the allele spends “in transit”, because while very rare, the allele will be mostly in heterozygotes, and so each generation provides another opportunity for a different haplotype to recombine closer to the transiting B allele. A large, linked haplotype may still arrive and fix in the new patch, in which case the haplotype has literally hitchhiked across geography! Fig 6 shows a simulation of such an event, including the lineage that founds an adaptive population on a second patch, and the length of the haplotype shared between this lineage and one in the original patch. The time that the lineage is outside the region where the B allele is common (dark red in Fig 6), the haplotype that accompanies it is broken down rapidly. After the lineage establishes on the patch, the rate of decay of the haplotype is slowed significantly, since most others with which it recombines have similar haplotypic backgrounds. Above we argued that a good model for this transit is a single Brownian trunk lineage surrounded by a cloud of close relatives of random size K whose chance of surviving until t generations is 1 − ke (t). Consider a single such family, and let τ be the (random) time at which it first hits the new patch, with τ = ∞ if it never does. We assume that the length of the hitchhiking segment depends only on the time spent by the trunk lineage of the first successful migrant “in transit” between the patches, i. e. , τ conditioned on τ < ∞. In so doing, we neglect recombination between B alleles (since they are at low density in transit), and the possibility that more than one successful migrant family is in transit at once (so that faster migrants would be more likely to have arrived first). Since each generation spent in transit provides an opportunity for recombination, if recombination is Poisson, the length of the haplotype (in Morgans) initially shared between populations on each side of the selected locus is exponentially distributed with mean τ. Therefore, if L is the length of hitchhiking segment on, say, the right of the selected locus, then P { L > ℓ } = E [ e - ℓ τ | τ < ∞ ]. (15) Computing this depends on 1 − ke (t), the probability a family survives until t. Since 1 - k e (t) ≃ e - s m t / E [ K ], we approximate the lifetime distribution of a family by an exponential with mean 1/sm. We can then use standard results on hitting times of d-dimensional Brownian motion that is killed at rate sm (see [45] 2. 2. 0. 1 and 4. 2. 0. 1). In particular, if the patch is circular with radius w and lies at distance R from the already adapted patch, then E [ e - ℓ τ ] = e - R 2 (s m + ℓ) / σ d = 1 K 0 ( (R + w) 2 (s m + ℓ) / σ) K 0 (w 2 (s m + ℓ) / σ) d = 2, (16) where K0 is a modified Bessel function of the second kind. We are interested in lineages that manage to reach the patch before being killed, i. e. , having τ < ∞, which occurs with probability P { τ < ∞ } = lim ℓ → 0 E [ exp (- ℓ τ) ]. To keep the expressions simple, in the remainder of this section we only compute quantities for one dimension. By Bayes’ rule, P { L > ℓ } = exp - R σ 2 (ℓ + s m) - 2 s m. (17) This can be differentiated to find that the expected transit time is E [ τ ∣ τ < ∞ ] = (R 2 s m / σ) × 1 / (2 s m) (18) and that Var [ τ ∣ τ < ∞ ] = (R 2 s m / σ) × 1 / (2 s m) 2. (A saddle point approximation provides an alternative route to these expressions.) The form of Eq (17) implies that if Y is an exponential random variable with rate R 2 / σ, then L has the same distribution as (Y + s m) 2 - s m. Furthermore, the expected length of shared hitchhiking haplotype is E [ L ] = σ 2 s m / R + σ 2 / R 2 (19) and Var [ L ] = 2 s m σ 2 / R 2 + 4 σ 3 2 s m / R 3 + 5 σ 4 / R 4. For two dimensions, asymptotics of Bessel functions show that L has the same tail behavior, but how other properties might change is unclear. As a rule of thumb, Eq (18) means that families who successfully establish move at speed 2 s m σ towards the new patch—if sm, the strength of the selection against them, is smaller, the need to move quickly is less imperative. Then, Eq (19) means that the haplotype length is roughly the length one would expect given the mean transit time, so more weakly deleterious transiting alleles arrive with them shorter haplotypes. Note that we have become used to seeing σ divided by 2 s m, rather than multiplied; it appears here because σ / 2 s m is a length scale, and to convert it to a speed this is divided by 1/2sm. Coat color in the rock pocket mouse Chaetodipus intermedius is a classic example of local adaptation [14,15]. These mice live on rocky outcrops throughout the southwestern United States and northern Mexico, and generally have a light pelage similar to the predominant rock color. However, in some regions these mice live on darker substrates (e. g. , old lava flows), and in these patches have much darker pigmentation, presumably to avoid visual predators. Some of the largest patches of dark rock range from 10km to 100km wide and lie 50–400km from each other, and dark-colored populations of C. intermedius have been found on many of these. (However, patches of all sizes occur across all scales in a heterogeneous manner.) [46] demonstrated that on one of these flows (Pinacate), much of the change in coloration is due to an allele at the MC1R locus. This dark allele differs from the light allele by four amino acid changes, and has a dominant or partially dominant effect depending on the measure of coat color. The Pinacate allele is not present in a number of other populations with dark pelage, suggesting these populations have adapted in parallel [17,46]. However, [47] reasoned that, elsewhere in the range, multiple close dark outcrops may share a dark phenotype whose genetic basis has been spread by migration despite intervening light habitat. A key parameter above was the dispersal distance divided by the square root of strength of selection against the focal allele between patches, σ / s m. [48] studied the frequency of the dark MC1R allele and coat color phenotypes, at sites across the (dark) Pinacate lava flow and at two nearby sites in light-colored rock. On the lava flow the dark MC1R allele is at 86% frequency, while at a site with light substrate 12km to the west (Tule), the frequency is 3%. The dark allele was entirely absent from Christmas pass, a site with light substrate 7km north of Tule, and 3km further from the lava flow. In the other direction, the dark MC1R allele was at 34% at a site with light substrate 10km to the east of the flow (O’Neill). Note that such apparent asymmetry is expected, since as noted above the migration–selection equilibrium can be highly stochastic. These numbers give us a sense of a plausible range of the parameters. Assuming the dark allele is at 50% frequency at the edge of the lava flow, we can fit Formula (3) to these frequencies (similar to [48]). Doing this for Tule we obtain σ / s m ≈ 3 km, and for O’Neill, σ / s m ≈ 30 km, giving us a range of cline widths. We also need estimate of sm. Using σ ≈ 1km [49,50], these cline widths imply that sm = 1/9 and 1/900 are reasonable values. The mutational target size μ for the trait is unclear. While the Pinacate dark haplotype differs from the light haplotype at four amino acid residues, it is likely that not all of these changes are needed for a population to begin to adapt. Also, there a number of genes besides MC1R at which adaptive changes affecting pigmentation have been identified in closely related species and more broadly across vertebrates [51]. To span a range of plausible values, we use a low mutation rate of μ = 10−8 (a single base pair), and a high mutation rate μ = 10−5 (a kilobase). Finally, we set A = 100km2 (roughly the size of the Pinacate patch). In Fig 7 we show the dependence of the probability of parallel mutation on the distance between lava flow patches using these parameters, showing that parallel mutation should become likely over a scale of tens to a few hundred kilometers between patches. Given the large selection coefficient associated with the dark allele on the dark substrate, we expect the initial haplotype associated with either a new mutation or migrant allele to be large. Fig 7 also shows how long the founding haplotype shared between populations is expected to be, from Eq (19). The initial length can be quite long between geographically close patches (tens of kilometers). However, for the wider cline width (σ / s m = 30 km), adaptation by migration can still be likely for patches 100km apart, but the shared basis may be hard to detect, as the length of shared haplotype can be quite short.
We have focused on relative rates of adaptation, since in applications where adaptation has occurred, the question is whether adaptations in distinct patches have appeared independently or not. However, any adaptation that does occur may have to make use of standing variation, if mutation rates are low. The case of a panmictic population was studied by [8], and we study the case of a continuous, spatial population in [54]. If parallelism in local adaptation of the sort we study here is due to standing variation rather than new mutation, then the dynamics of adaptation should not depend strongly on migration patterns (but the initial spatial distribution of standing variation may). We have mostly ignored the issue of dominance by dealing with essentially haploid models, and appealing to the fact that the dynamics we study occur where the mutation is rare, and hence mostly present only in heterozygotes. Our results should hold as a good approximation to dominant and partially dominant alleles (with sm the selection against heterozygotes). If, however, the mutation is recessive, then it is essentially neutral where rare, and so would encounter much less resistance to spreading between patches. The shape of the cline obtained is given by [24]. This is counteracted, however, by the increased difficulty with which the mutation would establish once arriving in a new patch, if the beneficial effect is also recessive. As such it is not clear what our intuition should be about the contribution of recessive alleles to adaptation via migration. Further work is needed to put empirical observations of local adaptation by recessive alleles in a theoretical context. It is similarly unclear how the model should extend to a polygenic trait. To provide context for the results on shared haplotype length in section Length of the hitchhiking haplotype it is important to also understand the process by which haplotypes are whittled down within patches. The initial haplotype that sweeps within a patch will be dispersed over time by recombination. Likewise, the haplotype that is shared between patches coadapted by migration will also break down (Eq 19). However, a long time after the initial sweep, we may still expect to find individuals within the patch sharing longer haplotypes around the selected locus than with individuals elsewhere, since selection against migrants decreases mean coalescence times within the patch near the selected locus. The literature on clines (e. g. , [44]) has important information, but more work is needed to provide robust estimates for these processes. Questions about the genomic length-scale of signals of sweeps shared by migration have also been addressed in discrete population settings [55,56], reviewed in [57]. This work has shown that the length of the shared swept haplotype is often significantly shorter than the sweep within each patch, resulting in a pattern of shoulders of elevated FST between adapted populations some distance away from the shared selected allele. It would be of interest to see how similar patterns can arise in a continuous population setting, as a way of uniting these results. We have also ignored the possibility of very long distance migration, instead focusing on local dispersal (hence Gaussian by the central limit theorem). However, dispersal distributions can be very heavy tailed, with a small fraction of individuals moving very long distances indeed [22,58]. In addition, over long time-scales, very rare chance events (mice carried off by hurricanes and the like; [59,60]) could play a role in spreading migrant alleles if adaptation by other means is sufficiently unlikely. Such tail events could greatly increase the probability of shared adaptation above that predicted by our model. Furthermore, if adaptive alleles do move between distant patches via rare, long distance migration then they will be associated with a much longer shared haplotype than predicted by Eq (19). As such, we view our results as a null model by which the contribution of long distribution migrants to adaptation could be empirically judged. We have studied circular patches of habitat at long distances from each other. Real habitat geometry can be much more complex, e. g. , with archipelagos of patches of varying sizes, or patches connected by long, skinny corridors, for instance. The work of [61] comes closest to a general theory of balanced polymorphisms in such habitats. It is possible that our techniques could be applied in their much more general setting, as both are based, fundamentally, on branching process approximations. It is also interesting to think about the probability of convergent adaptation to continuously varying environments, e. g. replicated environmental clines. The falling cost of population genomic sequencing means that we will soon have the opportunity to study the interplay of adaptation with geography and ecology across many populations within a species. Our work suggests that even quite geographically close populations may be forced to locally adapt by repeated, convergent, de novo mutation when migration is geographically limited and selective pressures are divergent. Thus, systems where populations have been repeatedly subject to strong local selection pressures may offer the opportunity to study highly replicated convergent adaptation within a similar genetic background [1]. Such empirical work will also strongly inform our understanding of the ability of gene flow to keep ecologically similar populations evolving in concert [62]. Our results suggest that adaptation to shared environments is certainly no guarantee of a shared genetic basis to adaptation, suggesting that rapid adaptation to a shared environment could potentially drive speciation if the alleles that spread in each population fail to work well together [63].
First we briefly describe the simulations we used for illustration and validation (the R code used is provided in S1 Scripts and at http: //github. com/petrelharp/spatial_selection). We simulated forward-time dynamics of the number of alleles of each type in a rectangular grid (either one- or two-dimensional) of demes with fixed size N. Each generation, each individual independently chose to reproduce or not with a probability r depending on her type and location in the grid; locally beneficial alleles were more likely to reproduce. Each extant individual then either remained in the same location with probability 1 − m or else migrated a random number of steps in a uniformly chosen cardinal direction; for most simulations m = 0. 2 and the probability of migrating k steps was proportional to 2−k for 1 ≤ k ≤ 5. In 2D, diagonal steps were also used. This gave us values of σ = 0. 95 deme spacings in 1D and σ = 0. 74 in 2D. Once migrants were distributed, each deme was uniformly resampled back down to N individuals. (Although we described the simulation in terms of individuals, we kept track only of total numbers in an equivalent way.) The base probability of reproduction in each generation in simulations for type b alleles was r = 0. 3; this was then multiplied by 1 + s to get the probability of reproduction for type B, where the value of s is either sm or sp depending on the individual’s location. This determines the values of sm and sp reported in the figures, and do not depend on the basic rate of reproduction. However, to obtain values for sp and sm when comparing theory to simulation, we computed the rate of intrinsic growth, i. e. , the s so that the numbers of B alleles when rare would change by est after t generations in the absence of migration. (The resulting values are close to the first notion of s, but give better agreement with theory, which uses the second definition.) To sample lineages, we first simulated the population dynamics forwards in time, then sampled lineages back through time by, in each generation, moving each lineage to a new deme with probability proportional to the reverse migration probability weighted by the number of B alleles in that deme in the previous time step. If more than one lineage was found in a deme with n alleles of type B, then each lineage picked a label uniformly from 1 … n, and those picking the same label coalesced. Since reproduction is Poisson, this correctly samples from the distribution of lineages given the population dynamics. When rare, copies of a new mutant allele are approximately independent and experience a uniform selective benefit; and can therefore be treated as a branching process. Furthermore, whether or not a new, beneficial mutation establishes or is lost to demographic stochasticity is determined by this initial phase where it is rare. Fortunately, the probability that a branching process dies out can be found as a fixed point of the generating function of the process [38]. Therefore, we calculated explicitly the generating function for a spatial branching process with nearest-neighbor migration on a one-dimensional lattice and a Poisson number of offspring with mean 1 + s, where s could vary by location, and iterated this forward to convergence to obtain 1 − p (x), the probability a single mutation appearing at x would fail to establish. We considered two situations: where s is a step function, and where it has a linear transition. These solutions are shown, and parameters described, in Fig 8. Here (and at other parameter choices) we see that the probability of establishment p (x) goes to the equilibrium value (approximately pe = 2s/ξ2) within the patch; the transition is fairly symmetrical about the edge of the patch, even if the edge of the patch is not sharp. Additional experimentation indicated that the fit remains equally good for other parameter values, even if migration can move further than one deme and offspring numbers are not Poisson. This lends credence to our approximation that the integral of p (x) over the entire range is close to pe multiplied by the area of the patch. For completeness, and clarity as to the scalings on the relevant parameters, here we provide a derivation of the differential equations referred to above Eq (3), and establish the asymptotics given in that equation. One route to the “equilibrium frequency” of the allele outside the range where it is advantageous is as follows; see [9] or [29] (or [64] or [65] or [24]) for other arguments in equivalent models, and see [66] and/or [67] for a general framework for the stochastic processes below. Suppose that the population is composed of a finite number of small demes of equal size N arranged in a regular grid, and that selection (for or against) the allele is given by the function s (x), with x denoting the spatial location. Each individual at location x reproduces at random, exponentially distributed intervals, producing a random number of offspring with distribution given by X who then all migrate to a new location chosen randomly from the distribution given by x + R, where they replace randomly chosen individuals. If x + R is outside of the range, then they perish. Each individual’s time until reproduction is exponentially distributed: the reproduction rate is 1 if it carries the original allele, or is 1 + s (x) if it carries the mutant allele. Suppose that the number of offspring X has mean μ; the variance of X will not enter into the formula (but assume X is well-behaved). Also suppose that the migration displacement R has mean zero and variance σ2; in more than one dimension, we mean that the components of the dispersal distance are uncorrelated and each have variance σ2. Let Φ t N (x) be the proportion of mutant alleles present at location x at time t, and Φt (x) the process obtained by taking N → ∞ (which we assume exists). Denote by δx a single unit at location x, so that e. g. Φ t N + δ x / N is the configuration after a mutant allele has been added to location x. For 0 ≤ ϕ ≤ 1, we also denote by X ¯ ϕ the random number of mutant alleles added if X new offspring carrying mutant alleles replace randomly chosen individuals in a deme where the mutant allele is at frequency ϕ (i. e. hypergeometric with parameters (X, Nϕ, N (1 − ϕ) ) ); similarly, X ˜ ϕ is the number lost if the new offspring do not carry the allele (i. e. hypergeometric with parameters (X, N (1 − ϕ), Nϕ) ). (We like to think of Φ t N as a measure, but it does not hurt to think of ΦN as a vector; we aren’t providing the rigorous justification here.) Then the above description implies that for any sufficiently nice function f (Φ) that ∂ ∂ t E f (Φ t N) = N ∑ x E (1 + s (x + R) ) Φ t N (x + R) f Φ t N + X ¯ Φ t (x + R) N δ x - f (Φ t N) + E 1 - Φ t N (x + R) f Φ t N - X ˜ Φ t (x + R) N δ x - f (Φ t N) (20) = μ ∑ x E ∂ ϕ (x) f (Φ t) Φ t (x + R) - Φ t (x) + s (x + R) Φ t (x + R) (1 - Φ t (x) ) + O 1 N. (21) In the final expectation, R and Φ are independent. This follows by taking first-order terms in 1/N in the Taylor series for f, and the fact that E [ X ¯ ϕ ] = ϕ μ and E [ X ˜ ϕ ] = (1 - ϕ) μ. We can see two things from this: First, since this is a first-order differential operator, the limiting stochastic process Φ obtained as N → ∞ is in fact deterministic (check by applying to f (Φ) = Φ (x) 2 to find the variance). Second, if we want to rescale space as well to get the usual differential equation, we need to choose Var[R] = σ2 and s (x) to be of the same, small, order; this is another way of seeing that σ / s is the relevant length scale (as noted by [9]). More concretely, suppose that the grid size is ϵ → 0, that Var[R] = (σϵ) 2, and that the strength of selection is s (x) ϵ, suppose that Φt (x) is deterministic and twice differentiable, and let ξ (t, x) = Φt/ϵ (x); then the previous equation with f (Φ) = Φ (x) converges to the familiar form: ∂ t ξ (t, x) = μ σ 2 2 ∑ k = 1 d ∂ x k 2 ξ (t, x) + s (x) ξ (t, x) (1 - ξ (t, x) ). (22) Here we have taken first the population size N → ∞ and then the grid size ϵ → 0; we could alternatively take both limits together, but not if ϵ goes to zero too much faster than N grows. One reason for this is that at finite N, the process Φt is an irreducible finite-state Markov chain with absorbing states at 0 and 1; therefore, the inevitable outcome is extinction of one type or another, which is not the regime we want to study. In one dimension, we are done (and discuss exact solutions in S1 Text); in higher dimensions, we are more interested in the mean frequency at a given distance r from a patch. If we take a radially symmetric patch centered at the origin (so s only depends on r), and let ξ (t, r) denote the mean occupation frequency at distance r, then the polar form of the Laplacian in d dimensions gives us that Eq (22) is ∂ t ξ (t, r) = μ σ 2 2 ∂ r 2 ξ (t, r) + σ 2 d - 1 2 r ∂ r ξ (t, r) + s (r) ξ (t, r) (1 - ξ (t, r) ). (23) A radially symmetric equilibrium frequency ξ (t, x) = q (∣x∣), with s (r) = −s < 0 for all r > r0, solves for r0 < r < ∞, ∂ r 2 q (r) + d - 1 r ∂ r q (r) - 2 s σ 2 q (r) (1 - q (r) ) = 0 (24) lim r → ∞ q (r) = lim r → ∞ ∂ r q (r) = 0 0 < q (r) < 1 (25) Since q (r) → 0 as r → ∞, so q (r) (1 − q (r) ) ≈ q (r), it can be shown that the true equilibrium frequency q is close, for large r, to the solution to ∂ r 2 u (r) + d - 1 r ∂ r u (r) - 2 s σ 2 u (r) = 0 (26) lim r → ∞ u (r) = lim r → ∞ ∂ r u (r) = 0. (27) This has general solution given by a modified Bessel function: using [68] 8. 494. 9, the general solution is u (r) = C ′ (r - r 1) (2 - d) / 2 K (2 - d) / 2 (r - r 1) 2 s / σ, (28) where C′ and r1 are chosen to match boundary conditions. Asymptotics of Bessel functions ([68], 8. 451. 6) then imply that q (r) ≈ C r (1 - d) / 2 exp - r 2 s / σ + O (1 / r), (29) where C is a different constant. Here we make a more precise argument to back up Expression (11) for the migration rate. The argument made above in section Heuristics applies to general migration mechanisms, since it relies only on a decomposition of the migrant families upon hitting the new patch; but it is also imprecise in subtle ways that are difficult to formalize. Here we take a somewhat different tack, supposing that it suffices to model the spatial movement of a migrant family by following only the motion of the “trunk” (i. e. , the red line in Fig 3), and supposing this motion is Brownian, with variance σ2 per generation. (Recall σ is the dispersal distance.) We then use facts about Brownian motion and branching processes, to compute more precise versions of Eqs (4) and (5). We are approximating the dynamics of the focal allele in the region further away than r0 from the patch as the sum of independent migrant families, each of whose dynamics are given by a spatial branching process (as depicted in Fig 3). Call B (r0) the region closer than r0 to the patch, and ∂B (r0) its boundary. Denote by γ (x) the mean rate of outflux of migrant families from a point x ∈ ∂B (r0), i. e. , the time-averaged density of individuals near a point x in ∂B (r0) that are the founders of new migrant families. Expressions (4) and (5) are a simple product of the “outflux of families”, when in fact they should be an integral of γ (x) over possible locations. However, it will turn out that the integrand of Eq (5) is well-approximated by a constant multiple of the integrand of Eq (4). First consider Eq (4) for the equilibrium frequency. Suppose that Z is one such spatial branching process as above in section The genealogy of migrant families, started at time 0 with a single individual at x, and write S for the new patch. The mean occupation measure of Z in the region S, which we denote u (x, S), can be thought of informally as the expected total number of offspring of a family beginning at x that ever live in S. Let q (S) = ∫S q (y) dy denote the total equilibrium frequency in S. This is decomposed in Expression (4) as the sum of mean occupation densities of a constant outflux of branching processes from ∂B (r0): ρ q (S) = ∫ ∂ B (r 0) γ (x) u (x, S) d x. (30) Now we will decompose u (x, S) under the assumption that the marginal distribution of the spatial motion of a single lineage is Brownian. Let Bt be a Brownian motion with variance σ2, and τ† an independent Exponential (sm) time. The mean occupation time of Z spent in a region S is, u (x, S) = E ∫ 0 ∞ Z t (S) d t (31) = ∫ 0 ∞ E [ Z t ] p t (x, S) d t (32) = ∫ 0 ∞ e - s m t p t (x, S) d t (33) = ∫ 0 ∞ P x { τ † > t & B t ∈ S } d t (34) = E x ∫ 0 τ † 1 S (B t) d t, (35) where P x gives probabilities for Brownian motion begun at x (i. e. , B0 = x), and likewise E x. If we define τS to be the hitting time of S by the Brownian motion B, and μS (x) to be the hitting distribution of ∂S by BτS conditioned on τS < τ†, by the strong Markov property this is equal to u (x, S) = P x { τ S < τ † } E μ S (x) ∫ 0 τ † 1 S (B t) d t (36) = P x { τ S < τ † } g (A), (37) where now E μ denotes expectations for Brownian motion for which the distribution of B0 is μ, and g (A) is defined to be the latter expectation, which does not depend on x if S is circular (with area A). This form we can now compare to the Expression (5) for the outflux of successful migrants. Consider the probability that a migrant family beginning with a single individual at x will ever establish in the new patch. It would be possible to analyze this probability directly, as in [29]; but here we take a simpler route, approximating this by the chance that the trunk hits the new patch, multiplied by the chance that at least one member of the family escapes demographic stochasticity and successfully establishes in the new patch. Write h (x, S) for the probability that the Brownian trunk hits the patch S before the family dies out, and f (S) for the chance that the family manages to establish in the new patch, given that it successfully arrives. (This is approximately independent of x.) As for q (S) above, Expression (5) is properly an integral against the outflux γ (x): λ mig (S) = ∫ ∂ B (r 0) γ (x) h (x, S) f (S) d x. (38) If we make the approximation that Z hits the new patch only if the trunk of Z does, and recall that 1 − ke (t) is the chance that Z survives for t generations, then h (x, S) ≈ ∫ 0 ∞ (1 - k e (t) ) P x { τ S ∈ d t } (39) = ∫ 0 ∞ e s m t (1 - k e (t) ) e - s m t P x { τ S ∈ d t } (40) = ∫ 0 ∞ e s m t (1 - k e (t) ) P x { τ S ∈ d t & t < τ † } (41) ≈ 1 E [ K ] ∫ 0 ∞ P x { τ S ∈ d t & t < τ † } (42) = 1 E [ K ] P x [ τ S < τ † ]. (43) Therefore, we have that h (x, S) E [ K ] ≈ P x [ τ S < τ † ] = u (x, S) / g (A). (44) Since the integrands in Expressions (30) and (38) only differ by a factor that (at least asymptotically) does not depend on the distance between x and S, we can obtain the migration rate by multiplying the equilibrium frequency by this factor. This factor will not depend on A, because although f (A) and g (A) in principle depend on the patch size (and geometry), the dependence is very weak. For instance, the width of a circular patch only very weakly affects the chance of establishment of a new migrant that appears on its edge, as long as the patch is large enough. By Eqs (30) and (38), ρ q (S) = ∫ ∂ B (r 0) γ (x) u (x, S) d x (45) ≈ g (A) E [ K ] ∫ ∂ B (r 0) γ (x) h (x, S) d x (46) = g (A) E [ K ] f (A) λ mig (S). (47) Once we show that g (A) ≈ 1/ (2sm), we will have arrived at the result, Eq (10). The function g (A) is the expected amount of time that a Brownian motion begun on the edge of a disk of area A is expected to spend inside the disk before τ†. This is integral of the Green function for the Bessel process of the appropriate order, so using [45], and letting w be the width of the patch, in d = 1, g (A) = ∫ 0 2 w / σ e - y 2 s m / σ 2 s m d y = (1 - e - 2 w 2 s m / σ) / (2 s m) (48) and in d = 2, by [68] 5. 56. 2, g (A) = ∫ 0 2 w / σ 2 y K 0 (y 2 s m) d y (49) = 1 s m ∫ 0 2 w 2 s m / σ y K 0 (y) d y (50) = 1 2 s m 1 - 2 w 2 s m σ K 1 (2 w 2 s m / σ), (51) where K0 and K1 are modified Bessel functions of the second kind. In either case, g (A) ≈ 1/ (2sm), which is the approximation we use in the main text. In the development above we need to approximate the integral of Eq (3) across the area occupied by the new patch, q (S) = ∫S q (x) dx. The precise answer depends on the shape and orientation of the patches; but we aim for a usable approximation, primarily in terms of the shortest distance from the old patch to the new patch, denoted R, and assuming R is large enough we can take Eq (3) as an equality. Since q (R) decreases as ∣x∣ increases, q (S) ≤ A q (R) (52) where A is the area of S. This will be a good approximation if S is small. In one dimension, if S has length ℓ, q (S) = ∫ R R + ℓ C e - u 2 s m / σ d u (53) = C e - R 2 s m / σ σ 2 s m 1 - e - ℓ 2 s m / σ (54) ≤ C e - R 2 s m / σ σ 2 s m. (55) In two dimensions, suppose that T is a rectangle enclosing S with one axis aligned towards the original patch and transverse width w. Then q (S) ≤ ∫ T q (x) d x (56) ≤ w ∫ R ∞ q (r) r d r (57) = w C ∫ R ∞ (r 2 s m σ) − 1 / 2 e − r 2 s m / σ d r (58) ≤ w σ 2 s m C π (R 2 s m σ) − 1 / 2 e − R 2 s m / σ d r. (59) If we absorb π into the constant C, this is of the form (width) × (σ / 2 s m) × q (R), i. e. , roughly the area, after replacing the length by σ / 2 s m. In both cases, the upper bound is q (S) ≤ A′ q (x), where A′ is the area of the parts of S that are no more than R + σ / 2 s m away from the original patch. Lower bounds could be obtained along similar lines.
Often, a large species range will include patches where the species differs because it has adapted to locally differing conditions. For instance, rock pocket mice are often found with a coat color that matches the rocks they live in, these color differences are controlled genetically, and mice that don’t match the local rock color are more likely to be eaten by predators. Sometimes, similar genetic changes have occurred independently in different patches, suggesting that there were few accessible ways to evolve the locally adaptive form. However, the genetic basis could also be shared if migrants carry the locally beneficial genotypes between nearby patches, despite being at a disadvantage between the patches. We use a mathematical model of random migration to determine how quickly adaptation is expected to occur through new mutation and through migration from other patches, and study in more detail what we would expect successful migrations between patches to look like. The results are useful for determining whether similar adaptations in different locations are likely to have the same genetic basis or not, and more generally in understanding how species adapt to patchy, heterogeneous landscapes.
Convergent Evolution During Local Adaptation to Patchy Landscapes