Genes as Cues of Relatedness and Social Evolution in Heterogeneous Environments

There are many situations where relatives interact while at the same time there is genetic polymorphism in traits influencing survival and reproduction. Examples include cheater-cooperator polymorphism and polymorphic microbial pathogens. Environmental heterogeneity, favoring different traits in nearby habitats, with dispersal between them, is one general reason to expect polymorphism. Currently, there is no formal framework of social evolution that encompasses genetic polymorphism. We develop such a framework, thus integrating theories of social evolution into the evolutionary ecology of heterogeneous environments. We allow for adaptively maintained genetic polymorphism by applying the concept of genetic cues. We analyze a model of social evolution in a two-habitat situation with limited dispersal between habitats, in which the average relatedness at the time of helping and other benefits of helping can differ between habitats. An important result from the analysis is that alleles at a polymorphic locus play the role of genetic cues, in the sense that the presence of a cue allele contains statistical information for an organism about its current environment, including information about relatedness. We show that epistatic modifiers of the cue polymorphism can evolve to make optimal use of the information in the genetic cue, in analogy with a Bayesian decision maker. Another important result is that the genetic linkage between a cue locus and modifier loci influences the evolutionary interest of modifiers, with tighter linkage leading to greater divergence between social traits induced by different cue alleles, and this can be understood in terms of genetic conflict.

This is equation (1) in the main text. To see the result, for N i ≥ 2, pick two random, distinct group members and note that the parent of the first has a probability of 1/N i of also being the parent of the other. The purpose of the two-generation process of group formation is to implement a simple means of producing variation in relatedness. Our general approach can readily handle other, similar demographic processes whereby a group is derived from a pool of dispersers. Social interactions. The offspring group members engage in a social interaction, for instance a public goods game, and produce offspring in proportion to the payoff in this game. An individual's phenotype z represents an investment (strategy) in the game, and we assume 0 ≤ z ≤ 1. The payoff to an individual with phenotype z in habitat i is a function of z and the average investmentz of the individual's group. For polymorphic populations the group compositions will vary. To study the invasion of a rare mutant, we need the expected payoff to a (randomly chosen) rare mutant player of the game with phenotype z , in a population where the resident phenotypes z 1 and z 2 occur with frequencies p i1 and p i2 (where p i1 + p i2 = 1). We write this as which is equation (2) in the main text. For convenience we will sometimes also use the notationw where the mutant z k is interpreted as a modification of z k , as well as which is the expected payoff in habitat i to the resident phenotype z k . The expectation in equation (S3) needs to be worked out for the particular social interaction in question, which is especially easy if w i (z,z) depends linearly onz. We make use of this in a convenient example of a public goods game: where W i is a baseline payoff, b iz is the reward (same for all group members) of the average investment in the group, and c i z 2 is the cost of the individual investment z, which is assumed to be quadratically increasing.
For a rare mutant z , founding groups that contain z will predominantly have only one mutant individual. Based on this, we can compute a conditional expectation of the group averagez as where the coefficient r i of relatedness is given by equation (S1). The expected payoff in equation (S3) is then We use the notation d ik for the derivative of this expected mutant payoff in habitat i with respect to z , evaluated at z = z k . From equation (S8) we get which is equation (3) in the main text. Sequence of events in the life cycle. For simplicity, we assume that individuals are haploid over most of the life cycle. However, to explore the evolutionary consequences genetic recombination, we introduce sexual reproduction by assuming that there is a brief sexual phase in the dispersal pool in a habitat, involving diploid individuals and crossing over to produce the haploid 'asexual' individuals that found the groups as described above. Mating is random with respect to the dispersal pool and occurs before the forming of groups in the habitat. In order to have a well-defined census point, we specify the population composition at a time after the sexual phase, when groups have been formed and the public goods game is about to start. The sequence of events in the life cycle, starting right after the census point, is as follows: (i) public goods game with offspring production, (ii) within-and between-habitat migration of the offspring, forming a dispersal pool in each habitat, (iii) mating and recombination, and (iv) the next episode of group formation, including one asexual generation. An entirely asexual life cycle is a special case where the rate of recombination is zero. This is effectively the same as removing the event (iii) in the life-cycle. Concerning migration, we assume that a proportion m ij of the offspring formed after playing the game in habitat j join the pool of dispersers in habitat i. We might have m 11 = m 22 = 1 − m and m 21 = m 12 = m, where m is a rate of between-habitat migration.
Genetic cues and population dynamics of a dimorphism. We study two (or more) loci, one being a genetic cue locus with alleles with values x k , and the other an epistatic modifier of the effects of alleles at the cue locus. We classify individuals according to their current habitat, i = 1 or i = 2, and their genetic cue, either x 1 or x 2 . An individual does not directly observe the habitat, but selects a phenotype based on the genetic cue. We can regard z as a function of x as a genotype-phenotype mapping. The resident strategy is to develop phenotype z k in response to the cue allele x k . For this dimorphism z 1 , z 2 , let n ik be the number of individuals in habitat i with phenotype z k , i.e., individuals possessing the cue allele x k at the census point. The number of individuals in habitat i is We assume that the number of groups, g i , that form in habitat i is constant, so the number n i = g i N i of individuals in a habitat at the census point is also constant. We also assume that the groups turn over in synchrony. The frequency p ik of z k in habitat i is then In situations where there are no mutant alleles at modifier loci, recombination between the cue locus and modifier loci does not change the distribution of phenotypes in the dispersal pool. We can then write the resident population dynamics as Here we used the notation φ i = n i m i1 (w 11 n 11 +w 12 n 12 ) + m i2 (w 21 n 21 +w 22 n 22 ) . (S13) We can interpret φ i as the expected number of founding group members deriving from a random individual in the pool of dispersers in habitat i (there are n i 'openings' in habitat i and m i1 (w 11 n 11 +w 12 n 12 ) + m i2 (w 21 n 21 +w 22 n 22 ) individuals compete for these 'openings'). Note also thatw jk and φ i in (S12, S13) depend on n ik (t), so the population dynamics is non-linear.
Because n i1 + n i2 = n i , we can use (S12) for k = 1 to describe the population dynamics of a cue dimorphism: This equation is suitable to iterate to find a population dynamical equilibrium of a cue dimorphism. Alternatively, we can use (S12) for k = 2 to describe the population dynamics of a cue dimorphism: This is equivalent to (S14), but can be convenient to examine the invasion of cue allele x 2 into a monomorphism of x 1 . Thus, we can use these dynamics to delineate the region in phenotype dimorphism space where z 1 and z 2 , corresponding to the cue alleles x 1 and x 2 , can coexist. We do this by examining when either of them can invade a monomorphism of the other. The corresponding z 1 , z 2 values then form the set of possible resident dimorphisms. Let λ 10 be the leading eigenvalue of the 2 × 2 matrix described by (S15) for n i2 = 0, and thus n i1 = n i (p i1 = 1, p i2 = 0). This matrix describes the first order dynamics of the invasion of x 2 into a monomorphism of x 1 , so x 2 can invade when λ 10 > 1. In the same way, x 1 can invade a monomorphism of x 2 when λ 01 > 1, where λ 01 is the leading eigenvalue of the 2 × 2 matrix described by (S14) for n i1 = 0, and thus n i2 = n i (p i1 = 0, p i2 = 1). The region of coexistence is thus given by min(λ 10 , λ 01 ) > 1.
For a numerical procedure, it is convenient to have analytic expressions of these leading eigenvalues, which are readily found for 2 × 2 matrices. Mutant modifier invasion. We are interested in the invasion of a mutant genotype-phenotype mapping, for which the genetic cue x k corresponds to the phenotype z k , into a resident dimorphism z 1 , z 2 . We focus on an epistatic modifier locus where a resident allele produces the phenotype z k when the allele x k is present at the cue locus. Let n ik denote the population dynamical equilibrium for such a resident dimorphism. For a rare mutant modifier allele, let n ik be the (small) number of mutants in habitat i that are linked to the cue allele x k (and thus have the phenotype z k ). We use (n 11 , n 21 , n 12 , n 22 ) as a vector of dynamic variables describing the mutant invasion into a resident z 1 , z 2 dimorphism, and we assume that n ik n i .
The population projection matrix K is 4-dimensional, and we write its elements as K ikjl , so that where Here, h ikl is the probability that a mutant modifier allele at the census point in habitat i is linked to cue allele x k , given that its 'ancestor' at the previous census point was linked to cue allele x l . The recombination frequency between the genetic cue and modifier loci is ρ. Considering the inheritance of modifiers, we then have for the inheritance factor h ikl in (S18). In equation (S19), p ik = m i1w1k n 1k + m i2w2k n 2k m i1 (w 11 n 11 +w 12 n 12 ) + m i2 (w 21 n 21 +w 22 n 22 ) is the frequency of individuals in the dispersal pool in habitat i (prior to mating) that have allele x k at the cue locus, and δ kl is the Kronecker delta. Because we assumed that the n ik correspond to a resident population dynamical equilibrium, we find using (S12) thatp ik = p ik = n ik /n i . We can then write (S19) as which is is equation (9) in the main text The population projection matrix can be written using four component 2-by-2 matrices: The components are given by Note that the dependence on z l only appears in K 1l and K 2l . For the case where the mutant is equal to the resident (z k = z k ), we use the notation Invasion fitness. Given a demographic equilibrium for a resident dimorphism z 1 , z 2 , produced by a resident modifier of the genetic cue locus dimorphism, the invasion fitness for a mutant modifier producing z 1 , z 2 is equal to the logarithm of the leading eigenvalue of the population projection matrix for the mutant invasion, which is equation (4) in the main text. Here, λ is the leading eigenvalue of the matrix K in (S17, S18). We know that F (z 1 , z 2 ; z 1 , z 2 ) = 0, i.e., the leading eigenvalue of the matrix K in (S24) is λ = 1. For an evolutionary equilibrium, we should have F (z 1 , z 2 ; z 1 , z 2 ) ≤ 0 for any mutant z 1 , z 2 .

Details of results
Reproductive value and selection gradient. For a demographic equilibrium, with population projection matrix K in (S17) with z l = z l , so that K = K in (S24), we have the leading eigenvector n = (n 11 , n 21 , n 12 , n 22 ) and a corresponding 'left eigenvector' v = (v 11 , v 21 , v 12 , v 22 ) of reproductive values. Let u be the eigenvector n normalized to sum to one, i.e., and normalize v such that the scalar product of v and u is equal to one. We can make use of the reproductive values in the interpretation of evolutionary equilibria by noting that we can express the derivative of invasion fitness with respect to z l , at z l = z l , as the matrix product This follows from standard perturbation theory of operators. Note that these derivatives with respect to z l correspond to derivatives with respect to the modifiers of the mapping from cue to phenotype. With notation from (S9), we can use (S18) to express the mutant derivatives of the elements of the population projection matrix as Direct fitness interpretation. The question is now if we can interpret polymorphic evolutionary equilibria as satisfying principles corresponding to the the direct fitness approach of social evolution theory. Conditional on 'observing' (being linked to) the cue allele x k that induces z k , a 'strategic decision maker' at a modifier locus has information about the habitat, and thus about the social environment. The prior probability to be in habitat i is The conditional probability of being in habitat i, given the observation of the cue x k , is then which is equation (6) in the main text. These conditional probabilities, given a genetic cue, provide the basis for an interpretation. From equation (S23) we see that the dependence on z l only occurs in the matrix components K 1l and K 2l . Let us use the notation v k = (v 1k , v 2k ) and u l = (u 1l , u 2l ). Equation (S27) can be written as Using the notation p l = (n 1l + n 2l )/(n 1 + n 2 ), we get that We can interpret this as the conditional expectation (given the 'observation' of the cue allele x l ) of the fitness effect of a small change in z l away from z l , measured in terms of reproductive value. For each observation of the cue, there are 4 possibilities of going through the population cycle, corresponding to to a combination of the outcomes of the migration and 'cue inheritance' events. Note that p l q jl = u jl and that p l has the interpretation of a 'dilution factor', taking into account the the phenotypic reaction to the cue allele x l is exposed to selection only in that part of the population that observes this cue. Let us use the notation V jl to denote the (expected) reproductive value of dispersing offspring of a player in habitat j that has observed the cue x l , i.e., which is equation (8) in the main text. We can now write (S33) as which is equation (7) in the main text. Invasion of modifier of monomorphism. We can apply the above to the special case of a monomorphism, i.e., z 1 = z 2 = z, in which case we only need to keep track of the habitat i, and not the cue allele x k , in the equations for the modifier invasion. Let n i be the (small) number of mutant modifier alleles in habitat i, which modify z to z . The population projection matrix is wherew i is defined in (S3). Let us also use the notationw i , as in (S5). For z = z we have and one readily finds that (n 1 , n 2 ) is a leading eigenvector of K 0 with eigenvalue λ = 1. Note that n i is constant in our model. We can normalize the eigenvector to (u 1 , u 2 ) with unit sum, so that u i = n i /(n 1 + n 2 ). The corresponding left where we introduced the notation R 12 for the ratio. Using the normalization v 1 u 1 + v 2 u 2 = 1 for the reproductive values, we have Just as for the dimorphism case, we can write the selection gradient as where v = (v 1 , v 2 ) and u = (u 1 , u 2 ). Introducing and using the vector of reproductive values from (S39), the selection gradient at the monomorphism z is then This expression can be used to compute monomorphic equilibria numerically.

Numerical approach.
The evolutionary equilibria shown in Fig. 2 and Fig. 3 in the main text were computed numerically using the selection gradient from equation (7) in the main text. For this computation, we developed a C++ program that follows a path of small steps through z 1 z 2 -space, each of which increase the invasion fitness in equation (4) in the main text, until reaching an accurate approximation of the equilibrium. The location of monomorphic equilibria, illustrated in Fig. 2 in the main text, was similarly determined using the selection gradient (S42). The coexistence regions in Fig. 3 in the main text were computed numerically using the approach given by equation (S16). The C++ program is available from the corresponding author upon request.

Individual-based simulations
In all simulations we used the genotype-phenotype mapping which is equation (5) in the main text. Here a 0 + a g x is a 'liability' to develop the phenotype z given the cue x, with parameters a 0 and a g . The genetic cue was always encoded by a single locus, but the parameters a 0 and a g had different genetic architecture in different simulations. Our general approach was to develop C++ programs that directly implemented the sequence of events in the life cycle, using pseudo-random numbers to handle stochastic events, such as recombination and mutation. Alleles were encoded as real numbers, restricted to a particular range for each locus. Mutation of an allele was implemented as a random increment to the allelic value, having a bi-exponential distribution with a particular standard deviation for each locus. In case a mutated value exceeded the range of a locus, the value was set to the corresponding limit of the range. For the results we report, simulations were run for 20000 full life cycles with a rate of mutation of 0.0010, to generate variation, followed by another 20000 full life cycles with a smaller rate of mutation of 0.0001. We used a total population size of 40000, which is fairly large, but there is still some random variation in average phenotypes from simulation to simulation. We illustrate this variation in figures by showing the mean ± SD over several replicate simulations. For the simulations reported in Fig. 2 and Fig. 4 A and B in the main text there were 3 loci: a genetic cue locus coding for x and one locus for each of a 0 and a g . The latter two loci were fully linked, and their rate of recombination with the cue locus was ρ. The range of allelic effects was [−1, 1] at the cue locus, [−2, 2] at the locus for a 0 and [0, 5] at the locus for a g . In Fig. 2 in the main text, evolutionary change at the locus for a 0 was prevented by setting the mutational increments to zero for this locus. The allelic value at the locus was initialized to a (fixed) value that would allow evolutionary adjustments in a g and a polymorphism in x to closely approximate the values obtained from setting the selection gradient, given by [7] in the main text, to zero. These fixed values were as follows. For ρ = 0 and m = 0.01, 0.05, 0.10 we used a 0 = −0.5185168, −0.3610248, −0.3082787, and for ρ = 0.5 and m = 0.01, 0.05, 0.10 we used a 0 = −0.5625004, −0.5375327, −0.4446206. The purpose of this approach, with suitable fixed values of a 0 , is to illustrate simulation outcomes where there is little scope for genetic conflict to produce polymorphism at modifier loci. Different fixed values of a 0 essentially correspond to different genetic constraints. On the other hand, for the simulations in Fig. 4 A and B in the main text, effects at the locus for a 0 were free to evolve, resulting in the same outcome in some cases, but also in the transfer of polymorphism to the locus for a 0 , as illustrated in Fig. 4 B in the main text.
For the simulation shown in Fig. 4 C in the main text, the genetic cue locus had effects with a range of [0, 1]. The parameter a 0 was the sum of positive effects at 5 loci, each with a range of [0, 0.4], and negative effects at another 5 loci, each with a range of [−0.4, 0]. The parameter a g was the sum of positive effects at 10 loci, each with a range of [0, 0.3]. All loci were unlinked.
For the simulation shown in Fig. 4 D in the main text, the genetic cue locus had effects with a range of [0, 1] and the parameter a g was the sum of positive effects at 10 loci, each with a range of [0, 0.3]. However, the parameter a 0 had a more complex genetic architecture. There was a sum of effects at 5 loci with a range of [0, 0.4], and a sum of effects at another 5 loci with a range of [0, 0.4], but these sums of effects were each in turn controlled by another locus that put a ceiling on the maximum total expression. So if ξ would be the the total expression without a ceiling ξ max , then the expression was η with the ceiling, where η is given as follows: η = ξ max 1 − exp(−ξ/ξ max ) .
This is meant to correspond to a situation where a regulatory locus limits the expression at other loci. The parameter a 0 was the difference between these two regulated sums of effects. Again, all loci were unlinked in the simulation. The C++ source code for the simulation programs is available from the corresponding author upon request.