Assumptions

We consider two populations: surface-dwelling and cave-dwelling. We are interested in determining when the cave population will evolve blindness, i.e. become mostly comprised of blind individuals, as has occurred in numerous natural systems. We first assume that the surface and cave populations do not experience drift (i.e. populations are of infinite size). Additionally, immigration from the surface population into the cave affects the allele frequency in the cave population, but emigration from the cave to the surface does not affect the surface population, as we assume that the surface population is significantly larger than the cave population. Generations are discrete and non-overlapping, and mating is random. We track a single biallelic locus, where B is the sightedness allele and where b is the blindness allele.

The frequency of b is denoted by Q∈ [ 0,1] in the surface population and q∈ [ 0,1] in the cave population. On the surface, we assume that blindness is strongly selected against, and Q is dictated by mutation-selection balance. These and all subsequent variables are described in Table 1.

Table 1 Terms and variables Full size table

Calculating the frequency of the blindness allele

Within the cave, the life cycle is as follows. (1) Embryos become juveniles and experience constant, directional selection with relative fitnesses of w bb =1+s, w Bb =1+hs, and w BB =1, where s≥0 and h∈ [ 0,1]. (2) Juveniles migrate into and out of the cave such that a fraction m of adults come from the surface and 1−m from the cave, where 0≤m≤1. (3) Adults generate gametes with one-way mutation, where 0≤u≤1 is the probability that a functional B allele becomes a non-functional b allele. (4) Gametes unite randomly to produce embryos. Given this life cycle, we calculate the allele frequency of the daughter generation (q ′) via standard equations:

$$\begin{array}{*{20}l} & q_{j}= \!\frac{(1+s)q^{2}+ (1+hs)q (1-q)}{(1+s)q^{2}+2(1+hs)q (1-q)+(1-q)^{2}} \quad \text{selection} \end{array} $$ (1a)

$$\begin{array}{*{20}l} & q_{a} = q_{j}(1-m) + Q m \qquad\qquad\qquad \text{immigration} \end{array} $$ (1b)

$$\begin{array}{*{20}l} & q^{\prime} = q_{a} + (1-q_{a}) u \qquad\qquad\qquad \;\;\; \text{mutation} \end{array} $$ (1c)

Analysis of the change in allele frequency

The change in allele frequency in one generation is Δ q=q ′−q. The first derivative of the dynamics is informative about the behavior of the model under the influence of the different parameters. Selection and mutation are directional forces, and increasing s or u increases Δ q for 0≤q≤1 (i.e. derivatives are non-negative). Increasing h causes selection to be more effective at low q, as rare b alleles are exposed to selection, but less effective at high q, as rare B alleles are sheltered from selection; increasing h increases Δ q if \(0 < q < {\left (1+\sqrt {1+s}\right)}^{-1}\) and decreases it if \({\left (1+\sqrt {1+s}\right)}^{-1} < q < 1\) (i.e. the derivative is positive below this threshold and negative above it). Migration harmonizes the allele frequency in the cave population towards the surface population allele frequency. Thus increasing m increases Δ q for low q and decreases Δ q for high q (i.e. the derivative is positive only when 0≤q<q z (h,s,Q)≤Q, where q z is a function describing a threshold). However, increasing Q increases Δ q for 0≤q≤1 (i.e. the derivative is non-negative).

Identifying equilibrium allele frequencies

The model we have developed is an example of migration-selection balance [26–28], extended to also include mutation. An equilibrium exists for this model when Δ q=0. For small s, there is only one equilibrium, and it is near 0. For large s, there is only one equilibrium, and it is near 1. Three equilibria will only exist for moderate levels of selection (Fig. 1). If s=m=u=0, all forces of evolution are eliminated and Δ q=0 for 0≤q≤1. A lower bound for any valid equilibrium is \(\frac {mQ(1-u)+u}{m(1-u)+u}\) (Proposition 1). An upper bound for any equilibrium is 1−m(1−u)(1−Q) (Proposition 2). Furthermore, it is important to note that

$$ Q \le \frac{mQ(1-u)+u}{m(1-u)+u} \implies Q \le \hat{q} $$ (2)

Fig. 1 As selection increases, the evolutionary dynamics of the cave population changes. When s is low (red line; s=0), there is only one equilibrium: near 0. As s increases (blue–brown lines, s=0.05,0.1,0.15,0.2, and 0.25) the local maximum (upper hump) increases and crosses the x-axis, producing three equilibria. When s gets high enough (pink line; s=0.3), the local minimum (lower valley) also crosses the x-axis, resulting in one equilibrium again. The location of the equilibria are marked using vertical lines at the bottom of the chart. For all curves m=0.01, h=0, u=10−6, and Q=0.01. The figure on the right is an enlarged view of a small part of the figure on the left Full size image

indicating that the equilibrium frequency in the cave population will be greater than or equal to the allele frequency in the surface population. Intuitively, this result is obvious as positive selection and one-way mutation only add to the frequency of the blindness allele in the cave.

Assuming s>0, the solution to Δ q=0 are the roots of the following cubic polynomial

$$ g(q) = A q^{3} + Bq^{2} + Cq + D = 0 $$ (3)

where

$$ \begin{aligned} A &= -s(1-2 h)\\ B &= s(1\,-\,m(1\,-\,u)(1\,-\,Q)\,-\,h(3\,+\,u\,-\,m(1-u)(1-2Q)))\\ C &= -(m (1-u)+u)+sh(1+u - m(1-u)(1-2Q))\\ D &= Q m (1-u)+u \end{aligned} $$

There are three possible roots of this equation, corresponding to three possible equilibria. Depending on the parameter values, Eq. 3 may have three real roots or one real root and two imaginary roots. While the values of the roots of this polynomial can be expressed analytically, these equations are too complex to be helpful for understanding the system. For simplicity, we will let \(\hat {q}\) represent any possible equilibrium, and \(\hat {q}_{a} \le \hat {q}_{b} \le \hat {q}_{c}\), stand for the roots of Eq. 3.

Protected polymorphism

Rather than tackling the equilibria directly, we first demonstrate that the cave population has a protected polymorphism. A protected polymorphism exists if the allele frequency moves away from both fixation and extinction, i.e. Δ q>0 when q=0 and Δ q<0 when q=1. For q=0, Δ q=Q m(1−u)+u and q=0 will be an equilibrium only if Qm=0 and u=0; otherwise Δ q>0 at q=0. For q=1,Δ q=−m(1−Q)(1−u) and q=1 will be an equilibrium if m=0, Q=1, or u=1; otherwise Δ q<0. Thus a protected polymorphism always exists except at the edge cases Qm=u=0,m=0,u=1, and Q=1. In biological terms, the cave population will be polymorphic despite directional selection for b if there is some immigration from the surface population and the surface population is polymorphic.

Validity of equilibria

An equilibrium is only valid in our model if it is real and between [ 0,1]; otherwise, it is not biologically interpretable in this system. Because there is a protected polymorphism, there will be either 1 valid, stable equilibrium, or 3 valid equilibria in a stable-unstable-stable configuration, depending on the parameter values. While we have not exhaustively determined the parameter ranges under which there will be only one valid equilibrium, we have determined that if h≥1/3 or if h<1/3 and \(s h > \frac {m(1-u)+u}{1 + u - m (1-u)(1-2Q)}\), there will be only one valid equilibrium (Proposition 3).

We can also estimate the amount of selection required such that g(q)=0 (Eq. 3):

$$ \begin{array}{l} s_{q}(m,h,u,Q) \\ \,\,\,=\frac{m (1-u) (q-Q)-(1-q) u} {q \left(q-q\left(q+m (1-Q) (1-u)\right) -h (1-q) \left(m (1-2Q) (1-u)-(1-2 q)-u\right)+q\right)} \end{array} $$ (4)

This equation is not valid for all m∈ [ 0,1]. If the migration rate is low, \(m < \frac {(1-q)u}{(q-Q)(1-u)}\), no level of selection will make q an equilibrium, as all equilibria will be greater than q. Similarly, if the migration rate is high,

$$m > \frac{(1-q) \left(h (1-2 q+u)+q\right)}{(1-u) \left(h (1-q) (1-2 Q)+q (1-Q)\right)} $$

no level of selection will make q an equilibrium, as all equilibria will be less than q.

Dynamics and the evolution of blindness

The dynamics of the evolution of the cave population depend on the parameter values and the starting allele frequency, q 0 . — Our model is likely well behaved, e.g. no limit cycles or chaotic behavior, even though we provide no formal proof of this. — If there is one equilibrium, then the frequency of b will evolve monotonically towards it, i.e. \(q_{t} \to \hat {q}\) as t→∞. If there are three equilibria, then the frequency of b will evolve monotonically to \(\hat {q}_{a}\) if \(q_{0} < \hat {q}_{b}\) and to \(\hat {q}_{c}\) if \(q_{0} > \hat {q}_{b}\).

When the cave population is founded, its initial allele frequency will likely match the equilibrium frequency on the surface (q 0 =Q). Because \(Q < \hat {q}\) (Eq. 2), the allele frequency in the cave population will increase due to selection until it reaches the lowest equilibrium, i.e. q ∞ = inf{q:0≤q≤1 and Δ q=0}. Whether blindness evolves in the cave population depends on whether q ∞ ≥q ∗, where q ∗ is a researcher-chosen threshold for determining that the cave population is a “blind” population. For example, q ∗=0.5 would specify that the blindness allele is the majority allele, and q ∗=0.99 would determine that the blindness allele is approximately fixed. We can also focus on phenotypes, and let a=q 2+2q(1−q)h measure the average blind phenotype in the cave population; then

$$a_{\infty} \ge a^{*} \implies q_{\infty} \ge \frac{\sqrt{a^{*}(1-2h)+h^{2}}-h}{1-2 h} $$

We define s ∗ as the minimum level of selection required for cave population to become blind, given the other parameters, i.e.

$$s^{*} = \inf \{ s : s > 0\ \text{and}\ q_{\infty} \ge q^{*} \gg Q\} $$

For simplicity, we will only consider values of q ∗ much higher than the surface allele frequency. If there is one equilibrium, \(\phantom {\dot {i}\!}s^{*} = s_{q^{*}}({m,h,u,Q})\); however, if there are three equilibria, q t will evolve to the lower equilibrium and q ∞ ≈Q≠q ∗ (typically). Thus selection needs to be strong enough such that there is only one equilibrium; therefore,

$$\begin{aligned} s^{*} \approx \inf &\left\{s : s > 0\ \text{and}\ s \ge s_{q^{*}}(m,h,u,Q)\right.\\ &\quad\text{and}\ \left.\Delta(s,m,h,u,Q) < 0 \right\} \end{aligned} $$

where Δ(s,m,h,u,Q) is the discriminant of Eq. 3. Figures 2, 3, and 4 plot analytical solutions for s ∗ based on different thresholds. When m≫u, the ratio s ∗/m is roughly constant such that if q ∞ ≥q ∗ then

$$ {\begin{aligned} \frac{s^{*}}{m} \ge \max &\left\{ \frac{q^{*}-Q}{q^{*}(1-q^{*})\left(q^{*}+h(1-2 q^{*})\right)},\ \frac{1-6Q}{h}\right.\\ &\;\left.+\frac{2 Q-2 \sqrt{Q^{2}+ h Q \left(1-3 h (1-3 Q)-6 Q\right) }}{h^{2}} \right\} \end{aligned}} $$ (5)

Fig. 2 The level of dominance of the blindness allele (h) affects the level of selection (s) required to produce blind populations. Each line represents how strong selection must be relative to migration (m) for blindness to evolve in the cave population for a given level of dominance (s ∗/m), where s ∗ is the minimum level of selection required for the cave population to become blind. Regions above the curves produce populations that are blind and regions below do not. Each panel contains a different condition for defining whether the cave population is blind. a For the blind allele to become the majority allele requires stronger selection when blindness is recessive (h=0) compared to when the allele for blindness is dominant. b For the blind phenotype to become the majority phenotype requires stronger selection when blindness is recessive compared to when the allele for blindness is dominant. c For the blind allele to become fixed requires stronger selection when the allele for blindness is dominant compared to when it is recessive. d For the blind phenotype to become fixed requires stronger selection when blindness is recessive. The curves were calculated analytically with u=10−6 and Q=0.01 Full size image

Fig. 3 The frequency of the blindness allele (Q) in the surface population affects the level of selection (s) required to produce blind populations. The format of this figure follows Fig. 2: a shows when the blind allele becomes the majority allele; b shows when the blind phenotype becomes the majority phenotype; c shows when the blind allele becomes fixed; d shows when the blind phenotype becomes fixed. As Q increases, the amount of selection required to evolve blindness in the cave population decreases. A surface population with a high Q can be considered “pre-adapted” to the cave. The curves were calculated analytically with u=10−6 and h=0.5 Full size image

Fig. 4 The mutation rate of the blindness allele (u) affects the level of selection (s) required to produce blind populations. The format of this figure follows Fig. 2. As u increases, the amount of selection required to evolve blindness in the cave population decreases for small m, and blindness will evolve regardless of selection. The curves were calculated analytically with Q=0.01 and h=0.5 Full size image

See Appendix for derivation.

The neutral-mutation hypothesis

If blindness evolves neutrally in the cave population (s=0), the equilibrium allele frequency will be governed by mutation-migration balance, i.e. \(\hat {q} = \frac {mQ(1-u)+u}{m(1-u)+u}\). Similar to s ∗, we can define a critical value m ∗, such that if m<m ∗, the cave population will evolve blindness.

$$m^{*} = \sup \{ m : q_{\infty} \ge q^{*} \gg Q\} = \frac{(1-q^{*})u}{(q^{*}-Q)(1-u)} $$

Clearly, if u=0, the cave population will not evolve blindness without the influence of selection (or genetic drift). However, a completely isolated cave (m=0) will evolve blindness if there is mutation (u>0). As the migration rate increases, the equilibrium allele frequency decreases such that if m>m ∗, the cave population will not evolve blindness. Similarly, increasing the mutation rate increases m ∗, allowing blindness to evolve for higher immigration rates and demonstrating the importance of mutation to the evolution of blindness when selection is weak.

Recessive blindness

If blindness is recessive (h=0), we can evaluate the dynamics of three equilibria in more detail. First, we will simplify our model by assuming that u≪1 such that 1−u≈1 and

$$ \Delta q \propto s q^{2}\left[1- q - m(1-Q)\right] + \left[Q m +u - q \left(m +u\right)\right] $$ (6)

Weak-selection approximation

If selection is weak, then an equilibrium exists near q=Q. We use a second-order Taylor series at q=0 to determine the upper bound on s for the presence of three equilibria (i.e. when selection is so strong that an equilibrium near Q does not exist). The second-order series allows us to determine the lower two equilibrium points, although this approximation is inaccurate as q increases. This approximation gives us

$$ \Delta q \approx s(1-m)q^{2} - (m+u)q + mQ + u $$ (7)

after assuming that 1−Q≈1. This equation has two roots, which are the lowest two of three total equilibria,

$$\begin{array}{*{20}l} \hat{q}_{a,1} &= \frac{m+u-\sqrt{(m+u)^{2}-4 s (1-m)(mQ+u)}}{2 s (1-m)}\\ \hat{q}_{b,1} &= \frac{m+u+\sqrt{(m+u)^{2}-4 s (1-m)(mQ+u)}}{2 s (1-m)} \end{array} $$

These two roots exist only if

$$ \begin{aligned} &0 < \sqrt{(m+u)^{2}-4 s (1-m)(mQ+u)} \implies\\ &s < \frac{(m+u)^{2}}{4(1-m)(mQ+u)} \end{aligned} $$ (8)

which provides us with an estimate of the upper bound on s for the presence of three equilibria.

The derivative of Eq. 7 is \(\frac {\text {d}\Delta {q}}{\text {d}q}\left (q\right) = 2 s (1-m) q - (m+u)\), and an equilibrium will be stable if \( -2 < \frac {\text {d}\Delta {q}}{\text {d}q}\left (\hat {q}\right) < 0 \). From this, it can be easily shown that \(\hat {q}_{a,1}\) is stable and that \(\hat {q}_{b,1}\) is unstable.

Strong-selection approximation

In order to determine the lower bound on s for the presence of three equilibria, we assume that selection is strong enough such that u/s≈0 and Q/s≈0. Therefore,

$$ \Delta q \propto -q\left[q^{2} - \left[1-m(1-Q)\right]q + m/s\right] $$ (9)

and the equilibria can be described as

$$ \begin{aligned} \hat{q}_{a,2} &= 0\\ \hat{q}_{b,2} &= \frac{1}{2}\left(1-m(1-Q)-\sqrt{\left[1-m(1-Q)\right]^{2}-\frac{4m}{s}}\right)\\ \hat{q}_{c,2} &= \frac{1}{2}\left(1-m(1-Q)+\sqrt{\left[1-m(1-Q)\right]^{2}-\frac{4m}{s}}\right) \end{aligned} $$

The latter two equilibria will exist only if

$$s > \frac{4m}{\left[1-m(1-Q)\right]^{2}} $$

which provides us an estimate of the lower bound for the presence of three equilibria.

The derivative of Eq. 9 is \(\frac {\text {d}\Delta {q}}{\text {d}q}\left (q\right) = -3 q^{2} + 2 [1-m\) (1−Q)]q−m/s, and it can be easily shown that \(\hat {q}_{b,2}\) is unstable and \(\hat {q}_{c,2}\) is stable.

Validity of approximations

By substituting \(\hat {q}_{a,1}\) and \(\hat {q}_{b,1}\) back into Eq. 6, we obtain \(\Delta q = - s \hat {q}^{2}\left (\hat {q}-Qm\right)\). Thus, Δ q≤0, which indicates that \(\hat {q}_{a,1}\) overestimates \(\hat {q}_{a}\) and that \(\hat {q}_{b,1}\) underestimates \(\hat {q}_{b}\). By substituting \(\hat {q}_{b,2}\) and \(\hat {q}_{c,2}\) back into Eq. 6, we find that \(\Delta q = Qm+u(1-\hat {q})\). Thus Δ q≥0, which indicates that \(\hat {q}_{b,2}\) overestimates \(\hat {q}_{b}\) and that \(\hat {q}_{c,2}\) underestimates \(\hat {q}_{c}\). However, the error in our approximations is slight (Fig. 5).

Fig. 5 Our recessive-blindness equilibria approximations are accurate. The approximations developed in this paper (solid lines) are a good fit for calculated values of selection (s) that result in equilibrium for a given frequency of the blindness allele (q; circles) using Eq. 3. The dashed lines are our approximate bounds for the existence of three equilibria (i.e. for small and large values of s, there is one equilibrium; for intermediate values of s there are three possible equilibria). Other parameters are m=0.01, u=10−6, and Q=0.01 Full size image

Dynamics

Based on these approximations, the dynamics of the recessive-blindness system can be summarized as follows. First, there are three possible equilibria: \(\hat {q}_{a} \approx \hat {q}_{a,1}\), \(\hat {q}_{b} \in \left [ \hat {q}_{b,1},\hat {q}_{b,2} \right ]\), and \(\hat {q}_{c} \approx \hat {q}_{c,2}\). Second, there are four possible equilibria configurations: 1, 2a, 2b, and 2c.

Case 1, \(\frac {(m+u)^{2}}{4(1-m)(mQ+u)} < \frac {4m}{\left [1-m(1-Q)\right ]^{2}}\): only one equilibrium exists, and it is stable. The population will always evolve towards it.

Case 2, \(\frac {4m}{\left [1-m(1-Q)\right ]^{2}} < \frac {(m+u)^{2}}{4(1-m)(mQ+u)}\): depending on the strength of s, this case may have one of three possible configurations:

Case 2a, \(0 \le s < \frac {4m}{\left [1-m(1-Q)\right ]^{2}}\): Only one equilibrium exists, \(\hat {q}_{a}\), and it is stable. The population will always evolve towards it.

Case 2b, \(\frac {4m}{\left [1-m(1-Q)\right ]^{2}} < s < \frac {(m+u)^{2}}{4(1-m)(mQ+u)}\): All three equilibria exist; \(\hat {q}_{a}\) and \(\hat {q}_{c}\) are stable, while \(\hat {q}_{b}\) is unstable. If the population starts below \(\hat {q}_{b}\), it will evolve towards \(\hat {q}_{a}\). If it starts above \(\hat {q}_{b}\), it will evolve towards \(\hat {q}_{c}\).

Case 2c, \(\frac {(m+u)^{2}}{4(1-m)(mQ+u)} < s\): only one equilibrium, \(\hat {q}_{c}\), exists, and it is stable. The population will always evolve towards it.

Furthermore if q 0 =Q, the selection-threshold for blindness to be established in the cave population is

$$ {\begin{aligned} s^{*} \approx \max \left\{ \frac{m(q^{*}-Q)-u(1-q^{*})}{q^{*2}\left(1-q^{*}-m(1-Q)\right)},\ \frac{(m+u)^{2}}{4(1-m)(mQ+u)} \right\} \end{aligned}} $$ (10)

where q ∗ is the allele-frequency threshold.

Additive blindness and multiple alleles

Next we investigate a model where blindness is due to many additive (h=0.5) loci of small effect. This model is motivated by the identification of 12 additive loci corresponding to the difference in eye phenotypes between cave and surface populations of A. mexicanus [14]. First, we will make the following assumptions: (1) there are k unlinked loci with two alleles (for sightedness and blindness), (2) in the cave population the fitness of an individual is \(1 + s \frac {x}{2k}\), where x is the number of blindness alleles the individual carries, and (3) m, u, Q, and q 0 are identical at each locus.

Because the forces of evolution are equivalent at every locus, they will evolve identically, and the change in allele frequency due to selection is

$$ q_{j} = q \frac{1+q s + (1-q) \frac{s}{2k}}{1+q s} $$ (11)

See the appendix for a derivation. Next, we will simplify our model by assuming that u≪1 such that 1−2k u≈1 and

$$ \begin{aligned} \Delta q \propto &-s(1 - (2k -1)m) q^{2} + (s (1 - m (1 - 2 k Q))\\ &- 2 k (m +u)) q +2 k (Q m + u) \end{aligned} $$ (12)

This has a single, stable, valid equilibrium:

$$\begin{aligned} \hat{q} = \frac{s (1 - m (1 - 2 k Q))- 2 k (m+u) + \sqrt{(2 k (m (Q s-1)-u)-m s+s)^{2}+8 k s ((2 k-1) m+1) (m Q+u)}}{2 s ((2 k-1) m+1)} \end{aligned} $$

which decreases a k increases. Furthermore, if m≫u

$$\frac{s^{*}}{m} \ge \frac{2 k (q^{*}-Q)}{q^{*}(1-q^{*})} $$

In summary, the frequency of blindness alleles will increase in the cave population until they reach equilibrium, and they will be majority alleles if s≥s ∗>4k m.

Finite-population simulations

Constant migration

Cavefish live in small populations and strong levels of drift may play a significant role in the evolution of blindness in cave species. To investigate the impact of drift on our recessive-blindness model, we simulated diploid populations of size N=1000 (based on population estimates by [3]) by modifying our life cycle (Eq. 1) to include a finite population:

$$\begin{array}{*{20}l} q_{j} &= \frac{(1+s)q^{2}+ q (1-q)}{(1+s)q^{2}+(1-q^{2})} & & \text{selection} \end{array} $$ (13a)

$$\begin{array}{*{20}l} q_{m} &= q_{j}(1-m) + Q m & & \text{immigration} \end{array} $$ (13b)

$$\begin{array}{*{20}l} q_{a} &\sim \text{Binomial}(q_{m}, 2N) / 2 N & & \text{drift} \end{array} $$ (13c)

$$\begin{array}{*{20}l} q^{\prime} &= q_{a} + (1-q_{a}) u & & \text{mutation} \end{array} $$ (13d)

Here the adult population consists of 2N alleles sampled with replacement from the post-immigration gene pool.

For every simulation, u=10−6 and Q=0.01. These values were chosen because they are believed to be reasonable estimates, and because we previously examined the impact of varying Q and u (Figs. 3 and 4). We further explain the impact of altering these choices in the discussion. We set q 0 =Q, varied s from 10−6 to 102, and varied m from 10−8 to 1.

We simulated 100 replicates for each combination of parameters; simulations were conducted for 10,000 or 5,000,000 generations. For each set of parameters, we recorded the average q ′ frequency across these 100 populations at specific time points.

Our simulation results for finite populations are qualitatively similar to our analytical results for infinite populations. For high migration rates, the average allele frequency is similar to the infinite model, except that drift allows some populations that have three equilibria to evolve blindness (Fig. 6 b). However, at low migration rates (Q m<u), populations have low average frequency of b at 10,000 generations, unless s>1. As immigration decreased, these populations became dependent on de novo mutations to produce b, which is a slow process. At 5 million generations, which is close to the estimated age of cavefish populations [32], the average allele frequency is a better match to the results from the the infinite-population model (Fig. 6 c); however, it differs in two respects. (1) When selection is ineffective (2N s<1), the average allele frequency reflects mutation-migration balance. And (2) when migration is low (4N m<1), the average allele frequency shows increased variation. Thus drift is the strongest force affecting the change in allele frequencies in the bottom left of Fig. 6 c. For N=100 (not shown), results are qualitatively similar, but stronger selection is required to overcome the stronger effects of drift present at smaller population sizes, which most often leads to loss of the rare blindness allele.

Fig. 6 Populations evolve blindness in the face of immigration only with the help of strong selection. a The equilibrium frequency of the blindness allele (q) for an infinite population, and b–d average frequencies of the allele after t generations in finite populations with either constant or episodic migration. For each combination of selection (s) and migration (m) we conducted 100 replicate simulations with fixed values of the mutation rate (u=10−6), frequency of the blindness allele in the surface population (Q=0.01), and q 0 =Q. Colors correspond to the frequency of the blindness allele for a given combination of s and m, where blue is high frequency (blindness evolved) and red is low (blindness did not evolve). The solid white line corresponds to the degree of selection required in the infinite population (a) to result in q ∞ >0.5 (\(s^{*}_{0.5}\)). The area between the solid and dashed lines corresponds to the region where three equilibria exist. If 2N s≪1, drift is stronger than selection, and if 4N m≪1, drift is stronger than migration. If m Q≪u, mutation is the primary force introducing copies of the blindness alleles to the cave population Full size image

Episodic migration

Because cave and surface populations may be connected intermittently due to flooding, we simulated periods of immigration followed by periods of isolation following a first-order Markov process. The probability of switching between from isolation to immigration or vice versa was 10% in each generation. Results for the intermittently connected simulations were nearly identical to previous simulations, with the exception that at high levels of migration and selection, drift was more effective in increasing average allele frequencies (Fig. 6 d).

Multiple loci

To determine the effects of drift with multiple loci, we implemented the following individual-based simulation:

$$\begin{array}{*{20}l} q_{j,i} &= \bar{w}^{-1} \sum_{a=1}^{N}{w_{a} q_{a,i}} & & \text{selection} \end{array} $$ (14a)

$$\begin{array}{*{20}l} q_{m,i} &= q_{j,i}(1-m) + Q m & & \text{immigration} \end{array} $$ (14b)

$$\begin{array}{*{20}l} q_{u,i} &= q_{m,i} + (1-q_{m,i}) u & & \text{mutation} \end{array} $$ (14c)

$$\begin{array}{*{20}l} q_{a,i}^{\prime} &\sim \text{Binomial}(q_{u,i}, 2) / 2 & & \text{drift} \end{array} $$ (14d)

where q a,i is the frequency of blindness allele at the i-th locus in the a-th individual, w a is the fitness of the a-th individual, and \(\bar {w}\) is the average fitness. Note that we use fecundity selection in this simulation to reduce its complexity. Simulations of 1000 individuals were run for 10,000 generations with u=10−6, Q=0.01, and a grid of s and m values. The number of loci was k∈{1,2,4,6,12}. For each parameter value, 100 simulations were run and several summary statistics were calculated: the average frequency of blindness alleles, the average fitness, the average phenotype, and the average genetic load in the cave population.

Our stochastic simulations agree with our deterministic results (Fig. 7 shows infinite and finite results for k=1 and k=12). As predicted by the deterministic model, increasing the number of loci increased the amount of selection required to evolve blindness in the cave population. This result is due to the fact that genes with smaller effect size show more genetic load due to migration of surface individuals into the cave. Even when migration was weak, a smaller effect size decreased the strength of selection relative to drift.