Appendix a Description of the Individual Based Simulations
Appendix A – Description of the individual based simulations
Each individual was given two allelic values of the ecological traitu(x1, x2) and two of the mating parameter m (m1, m2), corresponding to the two diploid loci. The phenotype (u, m) of each individual was simply the sum of the corresponding allelic values (u = x1 + x2, m = m1 + m2).
Mating: Each individual functioned as a hermaphrodite, i.e. both as a female and a male (this has no consequence for the evolutionary dynamics as long genetic drift can be ignored). Each female (i.e. each individual) chose a male randomly with weights according to her preference function, given by her m-value (eq. 15 in the main text). Each pair produced f offspring with genes according standard Mendelian inheritance.The loci recombined freely. With a given probability each allele was mutated by adding a normally distributed random number with zero mean and a given standard deviation (see figure legends for parameter values of different simulations).
Survival: The survival probability of each individual was given by its u-value and the fitness function (eq. 3 or 20, depending on the model) divided by f. The probability was truncated to the range [0, 1].
Appendix B – Calculations of the fitness gradients
It is assumed that the resident population has two symmetricalleles for the u trait (–x and +x) and is monomorphic for the mating trait (m). Thus there are three u phenotypes in the population, all with the same m-value, which I will label 1-3. Denote the number of genotype i individuals before mating Ni. Then the number of juveniles of each genotype after mating are given by:
, (i, j, k = 1..3)(B1)
where Q(i,j, k) is the probability that a mating between genotypes j and k will produce an offspring of type i,P(j, k) is the probability that a type j female will choose a type k male (eq. 15), and N is total adult population size. After survival, the number of adults of next generation of each genotype are given by:
,(B2)
where w is the appropriate fitness function (eq. 3 or 20), u is a vector of the u-values and nis a vector of the corresponding number of juveniles.
First of all, the equilibrium distribution of genotypes has to be found. In principle, equations (B1) and (B2) can be iterated until the population has reached an equilibrium. In practise, however, that can take a very long time. I used the Newton-Raphson method with a numerically calculated Jacobian matrix to speed up calculations. Iterations were stopped when the difference in numbers between two consequitive generations (eqs. B1 + B2) was smaller than 10-12 for all genotypes.
A rare mutant allele x' for the u trait will turn up in heterozygotes (x', –x) and (x', +x) with phenotypes: u1' = x' – x and u2' = x' + x. The recursion equation for the mutant genotypes is linear:
,(B3)
where the vector M(t) contains the number of adults of each genotype at time t and T is a matrix with elements Tij given by:
,(B4)
where Q'(i, j, k) is the probability that a mating between a mutant of type j and a resident type k will produce a mutant offspring of type i, P'(j, k) is the probability that a mutant female of type j will choose a resident male of type k, and P''(k, j) is the probability that a resident female of type k will choose a mutant male of type j. The equilibrium resident equilibrium genotype frequencies, N*, were taken from the numerical procedure described above.
Since the mutant recursion equation (B3) is linear, invasion fitness is given by the dominant eigenvalue of the transition matrix T. The fitness gradient was calculated by numerical differentiation: the invasion fitness of mutants x – x and x + x was calculated and the difference, divided by 2x, gives the fitness gradient. I used x = 10-6.
The fitness gradient of the mating character m was calculated in a similar fashion, only now the mutant occurs in three genotypes, one per u genotype. For the numerical differentiation, I used m = 10-4. The somewhat high value was necessary due to the flatness of the fitness landscape of the mating character.
Appendix C – allowing for asymmetric branching
The model with a continuous resource landscape
To justify the analysis in the main text and investigate the possible consequences of an asymmetric branching process, I calculated two new types of gradient landscapes. They can be thought of as alternative ways of slicing the three-dimensional trait-space spanned by the three parameters x1, x2 and m.Figs. 4-8 in the main text are slices along the diagonal surface x1 = –x2. Figure C1a below shows a gradient landscape of the two ecological alleles x1 and x2 keeping m fixed at zero (random mating). The thick arrows show principal directions of evolution and the grey thin lines are deterministic evolutionary trajectories, assuming the rate of evolutionary change is proportional to the fitness gradient times the number of allele copies present in the population (eq. 8). The trajectories show first of all that the diagonal is not particularly stable – trajectories tend to go in parallel to it instead of converging to it. Instead, any asymmetry developed early in the branching process seems to be adjusted at the end, when all trajectories converge to the convergent stable, symmetric, configuration (–x*, +x*) (Fig. C1a, open circle). However, it is also evident that all trajectories eventually end up in the same place. The symmetric configuration is the only convergent stable point and it is globally stable (save for extreme situations that lead to extinction). This conclusion was confirmed for a larger region of trait space (–2 < x < 2) than is shown here and with asymmetric mutation rates, such that mutations in one direction are 10 times as frequent as in the other direction. Further, the fitness gradient for assortative mating remains positive in a large part of the trait space shown here – it is positive for all combinations of (x1, x2) below the dotted line in the upper left corner.
Similar plots to the Figure below can depicted at increasing levels of mating assortment (m), but show very similar results. The exact position of the convergent stable point changes (cf. Fig. 4), but it remains globally stable. The boundary where mating fitness gradient changes direction moves down/right until it includes the convergent stable point (above m 1), as can be predicted from Fig. 6. At even higher levels of assortment the boundary moves back again, such that even stronger assortative mating is selected for at the convergent stable point (above m 3.1, cf. Fig. 6).
To complete the picture, I also produced gradient landscapes parallel to the diagonal surface x1 = –x2, but with a constant asymmetry: x1 = –x2 + . They looked surprisingly similar to the landscapes shown in Figs. 6 and 7, with only minor differences (not shown).
In conclusion, it seems likely that a population branching according to this model will drift away from the symmetric diagonal, but such deviations from symmetry will only yield minor changes to the overall ecological divergence as well as the evolution of mating assortment. Further, asymmetries will disappear towards the end of the branching process and the final outcome remains the same. It follows that the conclusions drawn from the analysis of this model along the symmetric diagonal are robust to at least moderate asymmetric deviations.
Figure C1. Gradient landscapes of the two models when the mating behaviour is fixed to random mating (m = 0), but allowing for asymmetric ecological alleles (x1 and x2). The thick arrows represent the principal directions of evolutionary change. Dashed line: the fitness gradient of x1 is zero. Dash-dotted line: the fitness gradient of x2 is zero. An increased mating assortment (m) is selected for below and to the right of the dotted line. Grey lines: deterministic trajectories. a) The model with a continuous, unimodal resource landscape. Parameter values: r = 1, K = 1, = 0.4 (same as Figs. 4-6). b) The model with two discrete resources. Parameter values: b0 = 0.6, K1,2 = 50000, r1,2 = 1, a0 = 0.0004, c = 0.1, a = 0.6 (same as Fig. 8a)
The model with two discrete resources
Fig. C1b shows a gradient landscape of the two alleles,x1 and x2, assuming random mating (m = 0). As before, the thick arrows show the principal directions of evolution and the grey, thin lines show deterministic trajectories. As in the previous model (Fig. C1a) there is a symmetric point where both gradients are zero (open circle). This point corresponds to full ecological divergence with two specialised homozygotes. However, it is only convergent stable if parameters are fixed to the diagonal x1 = –x2. If asymmetric configurations are allowed, a more complicated branching pattern emerges. Initially, close to the origin in the lower right corner, trajectories converge towards the symmetric diagonal, but as the ecological differentiation is more or less completed trajectories start to diverge from symmetry. This divergence from symmetry is initially due to a directional selection on heterozygotes away from a fitness minimum while the two homozygotes are still close to fitness maxima with little directional selection. As asymmetry builds up, the maladapted, most extreme, homozygote decreases in frequency and selection on the heterozygote becomes increasingly important for the most extreme allele. Eventually, an asymmetric solution is reached where the heterozygote is fully specialised on one resource and the least extreme homozygote is specialised on the other (closed circles).
During this process of branching to an asymmetric configuration, assortative mating is selected for only initially, inside the region marked with a dotted line (Fig. C1b). It seems as if the species is predestined to one of the asymmetric solutions and random mating. This conclusion is however contradicted by the simulation in Fig. 8a, which stays close to symmetry (not shown) and evolves to full speciation. To resolve this contradiction it is necessary to study gradient landscapes similar to Fig. C1b, but at above-zero levels of assortative mating (m > 0). An increased mating assortment decreases heterozygote frequency and reduces the importance of the directional selection on heterozygotes. For the parameter values chosen here the symmetric solution becomes locally convergent stable at m = 0.41 and the asymmetric solutions vanish at m = 0.64, which is still not even half way to full speciation (Fig. 8a). At m > 0.64 the gradient landscape for the two alleles is very similar to that of the previous model (Fig. C1a), with a globally convergent stable, symmetric configuration.