Second‐order cooperation: Cooperative offspring as a living public good arising from second‐order selection on non‐cooperative individuals

Switching rate between cooperating and non‐cooperating genotypes is a crucial social evolution factor, often neglected by game theory‐inspired theoretical and experimental frameworks. We show that the evolution of alleles increasing the mutation or phenotypic switching rates toward cooperation is in itself a social dilemma. Although cooperative offspring are often unlikely to reproduce, due to high cost of cooperation, they can be seen both as a living public good and a part of the extended parental phenotype. The competition between individuals that generate cooperators and ones that do not is often more relevant than the competition between cooperators and non‐cooperators. The dilemma of second‐order cooperation we describe relates directly to eusociality, but can be also interpreted as a division of labor or a soma‐germline distinction. The results of our simulations shine a new light on what Darwin had already termed a “special difficulty” of evolutionary theory and describe a novel type of cooperation dynamics.

Evolutionary research in the past decades was marked by a great interest in social evolution, and a generally richer and more nuanced understanding of the evolution of cooperation. All these advancements have been strongly driven and inspired by mathematical models and computer simulations of biological cooperation, in parallel with the emergence of socio-microbiology. Typical game theory-inspired simulations feature individuals with binary, cooperate versus not-cooperate phenotypes (Nowak and May 1992;Allen et al. 2012). When they include an evolutionary component, the mutation rate parameter typically specifies the fixed probability of switching from one behavior to another. In these simulations, a higher switching (mutation) rate is unsurprisingly less favorable to cooperation (Nowak and May 1992;Allen et al. [The copyright line for this article was changed on 31 December 2019 after original online publication.] 2012), which is easily explained by a lower probability of neighboring individuals having the same behavior (West et al. 2006(West et al. , 2007 and is consistent with experimental evidences (Harrison andBuckling 2005, 2007).
Here, we are interested in the case where the rate of switching between cooperators and non-cooperators can itself evolve. On one hand it seems reasonable to expect that a low cost of cooperation would select for a low switching rate from cooperators to non-cooperators due to second-order, indirect selection for cooperation, which has indeed been shown experimentally (Harrison and Buckling 2007). On the other hand, when non-cooperators are the dominant type in the population, the direction of selection on the loci controlling the switching rate toward cooperation is unclear. One can construct verbal arguments for selection in either direction, motivating this work. We primarily focus on qualifying the direction of selection at the locus controlling the switching rate from non-cooperators to cooperators in a regime where cooperators are counter-selected, and show that this selection creates a relevant social dilemma in and of itself.
The inability of classical models of the evolution of cooperation to represent the evolution of the mutation rate between cooperators and non-cooperators is part of a larger problem. The main modeling assumption is that genotype = phenotype = binary behavior. In the context of cooperation, many studies in microbes raise the issue of more complex encoding of cooperative behaviors. For example, the production and release of colicins (toxins produced by some E. coli strains to kill competitors) depends on a phenotypic switch (Cascales et al. 2007): only a small fraction of cells from a clonal population produce the colicins. These colicins provide a benefit (killing of competitors) to the isogenic non-producing cells, who are resistant to the colicin but do not produce it, at the expense of the producing cells, who commit lysis to release the colicins in the environment. Calling cooperators only the cells producing the colicin would not make sense, because non-producing cells are isogenic. Put differently, the cooperative trait whose evolution must be explained is the propensity to stochastically produce colicins, and the "cheaters" (non-cooperators) would be the strains who have a lower propensity for such stochastic switch. Genetically speaking, the cooperative allele is the one coding for this phenotypic switch, and thus any individual bearing this genetic system is a cooperator, independently of whether it will switch toward the colicin-production phenotype during its lifetime. A cheater is then an individual bearing a loss of function mutation in this system, often a deletion of the colicin production system, but that is still immune to colicins produced by other individuals and thus benefits from them.
Recently, Ackermann and collaborators (Ackermann et al. 2008;Diard et al. 2013) focus on Salmonella typhimurium, a pathogenic, gut-living bacteria with a similar extreme altruistic behavior. To exclude competitors, a small fraction of the population invades the mucosa and induces inflammation, triggering host's immune response and usually dying in the process. All other individuals reap the benefits of this cooperative behavior because the immune response kills most competitor species. This scenario raises the same conceptual issue: designating only individuals that express the secretion system as cooperators is too restrictive. The cooperative behavior is triggered by a phenotypic switch, which is genetically controlled and has evolved in the subpopulation that does not enter the mucosa.
Here we consider the general case of selection on the propensity of non-cooperators to stochastically generate cooperator offspring, either by a phenotypic switch-as in the two aforementioned examples-or by mutation. The central goal of our work is to identify and quantify the selection pressures that act on switching rate modifier alleles. To do so, we have developed an individual-based simulation system featuring the individuals playing prisoner's dilemma on a 2D grid, but including an evolvable mutation rate from cooperators to non-cooperators and from noncooperators to cooperators, that can as well be interpreted and implemented as a phenotypic switching rate. The ability to switch toward cooperative behavior by mutation may at first glance seem mechanistically very different than the ability to switch by regulation that we already discussed with the examples of colicins and virulence factors. However, these two phenomena are conceptually very similar and are subject to the same types of selection pressures. Selection on the propensity of a genotype to mutate toward a different phenotype is part of the broader field of second-order selection (Tenaillon et al. 2001;Woods et al. 2011). While such phenomena have not been much studied in the context of cooperation, there are multiple examples of selection acting on the ability to mutate toward a particular phenotype due to selection increasing the likeliness of a particular mutation (contingency loci) (Silverman 1979;Abraham et al. 1985;Moxon et al. 1994Moxon et al. , 2006Anjuwon-Foster and Tamayo 2017). Going one step further, Colizzi and Hogeweg (2014) raise the question of the ability of a single sequence to encode an ecosystem, where second-order selection would shape the mutational landscape such that different likely mutants could stably interact together. In summary, it has been recognized for a long time that natural selection can act on the probability of a genotype to mutate toward a sequence giving an alternative phenotype, in a similar way that it can act to produce regulatory networks permitting switching between different phenotypes without changing the genotype.
We show that the evolution of a high mutation (or phenotypic switching) rate toward cooperation when non-cooperative alleles dominate satisfies the usual requirements of a social dilemma. The crucial step in doing so is using a definition of fitness that does not only take into account the phenotypes of one individual and its neighbors, but also the expected payoff of the descendants they will generate. Namely, an allele causing a high switching rate toward a cooperative behavior is costly to the individual bearing it because it increases the likelihood that some of its offspring will suffer a cost for cooperation. At the same time, it is beneficial to others because it increases of the proportion of individuals performing the cooperative behavior in the neighborhood. This is a classical cooperation game, between two types of individuals, the first with a high switching rate and the second with a low one. To distinguish them from primary cooperators, we call the defecting individuals with high rate of switching to cooperating phenotype the second-order cooperators and think of the situation as the second-order cooperation dilemma. We found a large range of parameters where cooperative alleles do not invade, but secondorder cooperation alleles, ones that increase mutation rate toward cooperation, do.
The second-order cooperation dilemma, which we identify and explore here, is of the same nature as the classical, "firstorder" dilemma. However, our framework goes far beyond merely determining the quantitative differences between the two, parameter values, or even exact conditions under which some type of cooperation may be maintained. It highlights and establishes the importance of developing the social evolution theory at a different level, one that includes not just the behavior and phenotypic properties of the individual, but also of its offspring. Past work has shown that even in simple binary simulations, the apparent cooperative steady state is actually a complex equilibrium driven by the life cycle and dynamics of cooperative patches (Misevic et al. 2015). Mutation/switching rate is a key parameter determining such dynamics, and as such, it must be considered an integral part of any simulation system that aims at capturing the full range of the cooperative scenarios in nature. Moreover, as we show in our study, allowing the mutation rate to evolve leads to qualitatively different outcomes, extends the definition of the cooperating phenotype, and is crucial to understanding the ecology and evolution of social dilemmas.

Historical Narrative and Connections with Eusociality
Our work is primarily inspired by examples of cooperation found in the microbial world. However, precisely identifying the nature of the social, cooperative phenotype is a more fundamental and widespread problem. In this section, we draw broader connections between our simulations and the much greater subject of the evolution of eusociality. Eusociality was cited by Darwin as a complication in his theory (Darwin 1859), but he considered the existence of non-reproductive individuals easy to explain because of potential benefits at community level. Instead, the real problem Darwin saw in eusocial insects is frequently forgotten. He had trouble explaining the many phenotypic, including behavioral, differences that exist between the queen and the workers, since workers do not transmit their characteristics and thus could not have gradually acquired them (Herb 2014). In the words of Darwin (1859): How the workers have been rendered sterile is a difficulty; but not much greater than that of any other striking modification of structure [...]. But I must pass over this preliminary difficulty. The great difficulty lies in the working ants differing widely from both the males and the fertile females in structure, as in the shape of the thorax, and in being destitute of wings and sometimes of eyes, and in instinct [...]. If a working ant or other neuter insect had been an ordinary animal, I should have unhesitatingly assumed that all its characters had been slowly acquired through natural selection [...]. But with the working ant we have an insect differing greatly from its parents, yet absolutely sterile; so that it could never have transmitted successively acquired modifications of structure or instinct to its progeny. It may well be asked how it is possible to reconcile this case with the theory of natural selection? Today, this "special difficulty" of the theory of natural selection (Darwin 1859), one primarily arising from trying to explain insect eusociality, has an obvious and straightforward mechanistic explanation: epigenetic and phenotypic plasticity. In a way, finding the answer was difficult because the question that troubled Darwin and others since him was not formulated clearly. The inherited behavior, phenotype whose evolution must be explained, is not "infertility and propensity for brood care of the queen's eggs", but "differentiation of one part of the eggs into infertile individuals with a propensity for brood care of the queen's eggs"-similarly to the interpretation made in the case of colicinogenic E. coli or mucosa-invading S. typhimurium, where the inherited cooperative behavior is a certain rate of phenotypic switching. As explained by Darwin (1859): Thus I believe it has been with social insects: a slight modification of structure, or instinct, correlated with the sterile condition of certain members of the community, has been advantageous to the community: consequently the fertile males and females of the same community flourished, and transmitted to their fertile offspring a tendency to produce sterile members having the same modification.
More generally, and more recently, eusociality has been a serious point of friction between proponents of group and kin selection (Wilson and Hölldobler 2005;Foster et al. 2006;Nowak et al. 2010;Abbot et al. 2011;Liao et al. 2015). Part of the misunderstanding can precisely be traced to inclusive fitness theory being applied to different objects (Fletcher et al. 2006). Instead of trying to understand the evolutionary interest of the workers, or the workers' genes, in not reproducing, one could focus on trying to understand the evolutionary interest of the queen, or the queen's genes, in having a large proportion of offspring being sterile. In the context of eusociality, this idea is often referred as "maternal control" (Lehmann et al. 2008), and is very similar to our idea of studying the fate of an allele increasing the switching rate from non-cooperative to cooperative behavior. Since here we implemented a prisoner's dilemma where cooperation modifies fitness, our framework can represent cooperative behavior that are less extreme than eusociality, both qualitatively and quantitatively. Namely, the fraction of cooperative offspring can be smaller compared to the one observed in eusocial species, and the cooperative behavior we model does not necessarily imply suicide or total loss of reproduction. However, we show that these assumptions can be relaxed without causing a qualitative difference in results. The problem of inheriting a conditional or probabilistic cooperative behavior is much broader than eusociality, further justifying our simple and general simulation system.

Methods
We simulate the evolution of individuals playing prisoner dilemma within their 3x3 neighborhood on a square 100x100 lattice with periodic boundaries (toroidal). The genotype of each individual is represented by three binary loci that control: (1) cooperative behavior, C (cooperative) or D (non-cooperative), (2) probability of mutation from D to C, L DC or H DC for a low or high mutation rate, and (3) probability of mutation from C to D, L C D or H C D . The fitness of an individual is determined by playing a prisoner dilemma with each individual in its 3x3 neighborhood: where N i is the classical 3x3 Moore neighborhood of the individual i that includes the individual i itself, C k is 1 if individual k bears allele C and 0 otherwise. Only relative fitness to the neighborhood matters, so we reduce the parameter space to one dimension by fixing base f itness = c − 1/8 and bene f it = 9/8. In our simulations we will vary the last remaining free parameter in the fitness calculation, the direct cost of performing the cooperative behavior, c.
The populations were initialized by randomly picking, for each individual, one of the two possible alleles for each locus. The simulations are synchronous and at each time point, for each location on the lattice, we compete the focal individual with its eight neighbors to determine which of them is going to leave the offspring in that location, in the next generation. This means that selection does not act on lifespan (each individual is replaced at each generation) but on fecundity: the individuals from the neighborhood are competing to reproduce, leaving a descendant to replace the individual previously occupying a given site on the lattice.
The cooperative behavior is heritable (except for Fig. S3) and can change due to mutations, which happen during reproduction. The probability of the mutation in the first locus is controlled by the other two loci, while the probability of mutations for the second and third locus is extrinsically fixed and unchanging. See Figure 4 for a graphical representation of all the possible genotypes on the mutational landscape. When the H DC and H C D mutation rates take high values, our system may capture the behavior of a cooperative trait controlled by a phenotypic switch with epigenetic inheritance. We also implemented simulations where cooperative behavior is not heritable (Fig. S3). Further simulation details are given in the "Simulation setup" section of the supplementary materials (Text S8).
Our simulations are stochastic, both due to random mutations described above and the probability-based reproduction. Specifically, the probability of reproduction of an individual i with fitness f i depends on its relative fitness, normalized by the fitness of all the individuals in the neighborhood: where m is the strength of selection. When m = 0, there is no selection, every individual has equal probability of reproduction, and the population evolves only by genetic drift. When m = 1, selection is linear and the probability of reproduction is proportionate to fitness. When m = ∞, the fittest individual is always the one picked for reproduction, which corresponds to a strong deterministic selection. The relevance of this framework is further discussed in the "Choice and relevance of our modeling framework" section of the supplementary materials (Text S10).
Any individual with allele C on the first locus, meaning that its genotype could be represented as (C, * , * ), is further referred as a pure cooperators. In contrast, second-order cooperators are individuals that do not express the cooperation gene but have a significant number of pure cooperator offspring and have genotypes (D, H DC , * ).

COOPERATION
In a first set of simulations, we varied the two main parameters influencing cooperation, cost of being a cooperator, c, and intensity of selection, m. c directly impacts the fitness difference between cooperators and non-cooperators, and is thus the main parameter influencing the outcome of a social dilemma in a spatially structured population. m impacts the strength of the competition that decides what individual from the 3x3 neighborhood is going to divide and replace the focal individual. When m is small, the competition has a higher degree of stochasticity (in the extreme case m = 0, all individuals of the 3x3 neighborhood have the same probability of reproduction and the population is only subject to genetic drift). When m is big, the competition is more deterministic (in the extreme case m = ∞, only the fittest individual of the 3x3 neighborhood reproduces).
After 2000 generations of evolution, we can identify the region of the parameter space in which individuals of type (C, * , * ), the pure, first-order cooperators (blue bars, Fig. 1), are abundant, or more precisely, are present at a higher level than predicted by the mutation-selection balance values observed in scenarios that select against cooperation. As expected, the higher the cost, the lower the proportion of pure cooperators. The effect of the intensity of selection is less straightforward, but we see that high intensity of selection can compensate for the high cost and select for some cooperation even when the cooperation is moderately costly (e.g., for c = 0.4 and m = 50).

COOPERATORS ARE ABUNDANT
After focusing on the evolutionary fate of the first locus, the one controlling the cooperative behavior of the individual, here we analyze the evolution of the other two loci. We start by considering the region of the parameter space in which allele C is abundant (panels with blue background in Fig. 1). Two main patterns emerge: (1) the second locus, controlling the rate of switching from D to C, has about equal proportions of the two alleles (green bars, blue shaded panels, Fig. 1) and (2) the third locus, control-ling the rate of switching from C to D, is dominated by the L C D allele, that is low mutation rate from C to D (red bars, blue shaded panels, Fig. 1). Given that cooperation is favored and defection is rare, selection on the second locus that affects D allele is weak, which explains the observed equal proportions of both alleles. The low proportion of alleles for high switching rate from C to D, H C D , is consistent with the idea that a higher mutation rate impedes cooperation by making more non-cooperators (cheaters) appear at the middle of cooperator patches. As cooperation is favored for this set of parameters, the high mutation rate that would hinder cooperation is in turn selected against. We confirm the role of mutation rate modifier alleles by running two sets of control simulations where we enforce respectively low and high mutation rate. Namely, we allow the first locus to evolve, but fix the genotypes on the two other loci, so that all the genotypes are either ( * , L DC , L C D ) in the first or ( * , H DC , H C D ) in the second set of control simulations. We combine the results from both control simulations and plot the average final proportion of pure cooperators as bars in Figure 2. Purple bars represent the control simulations with the constant low mutation rate, cyan bars ones with the constant high mutation rate, and blue bars are the same as in the Figure 1, corresponding to simulations in which mutation rates are free to evolve. For almost all simulations in which cooperation is abundant, it remains so in both control simulations. When mutation rate is constant and low there is only a slight increase in the proportion of cooperators (comparison between blue and purple bars, panels shaded blue, Fig. 2), which is not surprising, given that the evolved mutation rate is typically low as well. However, when high mutation rate is enforced, the number of cooperators is decreased (comparison between blue and cyan bars, panels shaded blue, Fig. 2). These simulations confirm that a lower mutation rate is favorable to 'pure' cooperation (Harrison and Buckling 2005;Allen et al. 2012) and can be selected for this reason (Harrison and Buckling 2007).
Finally, we should note that there are parameters on the edge of the space in which C is abundant (for example c = 0.4, m = 50), for which the effect of mutation rate (difference between number of cooperators in "low mutation rate" control, purple bars, and in "high mutation rate" control, cyan bars) becomes high. There may even exist a small parameter range where mutation rate from cooperators to non-cooperators would be decisive not only for the second-order cooperation but the pure cooperation as well. The indirect benefit of cooperation would be nearly completely balanced out by its direct cost and higher mutation rate from cooperators to non-cooperators would be enough to tilt the scale and enable a full invasion of the population by non-cooperator (cheater) mutants.

COOPERATORS ARE NOT ABUNDANT
Here, we analyze the majority of the parameter space, characterized by the cooperative allele C being rare and the allele D being dominant at first locus. Interestingly, in a significant part of this range, the allele H DC has also fixed, indicating a high mutation rate from non-cooperator to cooperator type (green bars, panels shaded green, Fig. 1). Put differently, the dominant genotype is (D, H DC , * ) and while the first-order, constitutive cooperation is outcompeted, majority of the individuals are second-order cooperators -they do not cooperate themselves but their offspring have a high probability of doing so. It is easy to show that secondorder cooperation is subject to the usual cooperation dilemma. The cost of second-order cooperation is having a fraction of less fit (or effectively sterile) offspring, and the benefit for the neighbors is the higher chance of encountering (C, * , * ) individuals, who are the offspring of second-order cooperators. Indeed, in this part of the parameter space, the aforementioned control simulations with fixed mutation rates confirm that the number of pure cooperators (allele C) present in the population is higher when enforcing a high mutation rate than when enforcing a low mutation rate

drawn as thin black lines (solid or dashed) when the mutation rate is low (0.001) and thicker arrows when it is high, due to mutation rate modifier allele H . More specifically, green arrows indicate the higher mutation rate from non-cooperator to cooperator due to H DC allele (second-order cooperation) and red arrows indicate the higher mutation rate from cooperator to non-cooperator
due to H C D allele. For each genotype, the one allele that does not directly affect the genotype is written in gray. For example, for the genotype (C , H DC , H C D ), top back right vertex in the figure, the second allele has no effect since it specifies the mutation from D to C allele and the individual is already C .
Supporting our claim that in this parameter zone pure cooperators (C, * , * ) are mainly present as offspring of second-order cooperators (D, H DC , * ) and not as descendants of other pure cooperators, we performed additional simulations in which pure cooperators are sterile. The outcome of these simulations is presented on Figure S6. Second-order cooperation remains prevalent in the same parameter zone, and additionally in the parameter zone where pure cooperators used to be dominant. These simulations also confirm that our framework can represent biological scenario where the act of cooperation is lethal, such as release of colicins through cell lysis.
While we do not track lineages in our simulations, classical population genetics can also provide a more quantitative reasoning about the expected frequency of pure cooperators under a mutation-selection balance scenario where they are mostly not reproducing but present by mutations from second-order cooperators. Here we are not in a textbook scenario due to spatial structure and local competition, but we can still get an approxi-mate expectation. In the simulations represented on Figure 1, the high mutation rate toward cooperation conferred by allele H DC is 0.01. In a population where second-order cooperation is dominant (green-shaded parameter zone), classical population genetics predicts the equilibrium frequency of pure cooperators (C, * , * ) to be equal to the frequency of second-order cooperators multiplied by 0.01 (mutation rate) divided by the selection coefficient. Here pure cooperation is costly and quickly goes extinct, so the selection coefficient is close to 1. Because second-order cooperation is prevalent in this parameter zone, the frequency of second-order cooperators is close to 1. We thus expect the frequency of pure cooperators to be close to 0.01. These frequencies are hard to visualize on Figure 1 that uses a linear scale, but can be read from Figure 2 that uses a logarithmic scale. As expected, in the green shaded parameter zone the equilibrium frequency of pure cooperators is very close to 0.01, consistent with our hypothesis that these individuals are mainly descendants of second-order cooperators.
The dilemma created by the H DC allele can be analyzed more mathematically by extending the definition of fitness to take into account not only the direct reproductive success (number of offspring) of the focal individual, but also the reproductive success of the offspring themselves; as done in the "Extended fitness measure" section of the supplementary materials (Text S9).
To investigate the effect of spatial structure on second-order cooperation, we calculate the regression coefficient between alleles H DC /L DC and allele C in the baseline simulations presented in Figure 1 in order to quantify the spatial association between genotypes. The coefficient is simply the average number of cooperative neighbors that individuals of type (D, H DC , * ) / (D, L DC , * ) have, averaged between generations 1000 and 2000. We performed paired one-sided Student's t-tests to compare the two different assortments and found that the assortment between (D, H DC , * ) and C is stronger than the assortment between (D, L DC , * ) and C, whenever D is dominant (purple and green regions, Fig. 3). Conversely, when C is abundant, the two aforementioned assortments are not significantly different. The lack of difference is due to cooperators forming patches and reproducing, instead of being only present as descendants of (D, H DC , * ) individuals.
The immediate consequence of such spatial clustering in regions where D is dominant is that H DC individuals benefit on average more from cooperation than L DC individuals. The benefit for second-order cooperators coming from this specific assortment is equivalent to the benefit of spatial clustering for the pure cooperators in classical, one locus simulations of cooperation. In our setup, the pure cooperators can be seen as helper individuals analogous to a living public good or a soma, created by (D, H DC , * ) individuals. Intuitively, the H DC allele may decrease F1 fitness but increase F2 fitness: one part of first-generation descendants are doomed, but provide a benefit that is more likely to be directed toward other first-generation descendants (due to spatial structure), potentially increasing the total number of second-generation descendants.
We should note that in a large part of the parameter space, one of the mutation-rate modifier loci is mainly evolving by genetic drift. When cooperators are dominant, the non-cooperators are too few for L DC /H DC to have an impact (blue region). Similarly, when non-cooperators invade (panels shaded green or purple), the cooperators are too few for L C D /H C D to have an impact. This outcome is visually represented in Figure 4: for any dominant genotype, gray shaded allele marks the locus that is not under selection, but subject primarily to drift dynamics.
To confirm that a high mutation rate from D to C is a cooperative behavior and evolves via group/kin selection, we performed additional simulations with randomized neighborhoods: for each focal individual, instead of calculating its fitness based on the benefit it gets from the cooperative behavior of the individuals in its neighborhood, we calculate it based on nine individuals picked randomly in the whole population. This is not exactly equivalent to a well-mixed population in which selection would be defined globally, but this allows to have the benefit of the cooperative behavior randomly spread in the entire population and not directed toward the neighbors. In these randomized conditions, for all parameters, the (D, L DC , * ) individuals invade the population. As before, there was no selection on the locus controlling mutation rate from C to D, simply because there were not enough cooperative individuals, so the two alleles on the second locus are present at equal frequencies. These results are represented on Figure S1 for high mutation rate 0.01. They are visually straightforward: without local interactions, both first and second-order cooperators are outcompeted, and the population is dominated by pure non-cooperators (D, L DC , * ). This is consistent with our interpretation of second-order cooperation, making H DC a cooperative allele even though the individuals bearing it do not necessarily express cooperation themselves.
As an additional evidence that the allele H DC creates a social dilemma in D individuals, we analyze the ecological dynamics of a rare H DC mutant in a L DC population, and of a rare L DC mutant in a H DC population, in control simulations presented on Figure S2. For several parameter sets, we run at least 100 replicate simulations of invasion of a single (D, H DC , * ) (or (D, L DC , * )) individual in a population of 399 (D, L DC , * ) (or (D, H DC , * )) individuals, for 200 generations. These simulations were performed without mutations on L/H DC and L/H CD loci. For each individual in the initial population, the allele at L/H CD locus is chosen randomly. The invasion success measure is the average, among replicate simulations, of the maximal frequency reached by the rare allele. Within each panel, we plot this invasion success measure both for structured populations (similar to Fig. 1, S4, S5) and when randomizing neighborhoods (similar to Fig. S1). High mutation rate is 0.01 for top row, 0.1 for middle row, and 0.5 for bottom row, corresponding respectively to the simulations presented on Figures 1, S4, and S5. We also perform, for each parameter set, at least 100 simulations of the same invasion success measure for a neutral allele, allowing us to compare our results with what would be expected for an allele only subject to genetic drift. For each condition we performed a Welch's t-test to determine whether there is a significant difference between the invasion success of a rare L DC and a rare H DC . The results are visually straightforward: in a well-mixed population (randomized neighborhoods), the L DC allele always invades more than a neutral allele, which in turn always invades more than the H DC allele. In a structured population, depending on the parameters, H DC may sometimes invade, at a low -but higher than genetic drift -frequency.

SWITCHING
So far, the mutation rate in our simulation has varied between 0.001 (L alleles) and 0.01 (H alleles). These values may be appropriate if we consider that the switching between cooperation and defection phenotypes is driven by genetic mutations, but could be much higher if phenotypic switches are involved (Dubnau and Losick 2006). While the reasoning about second-order cooperation remains exactly the same, the higher switching rates may be relevant in a broader range of natural scenarios. Using the same setup we test for the presence and significance of the second-order cooperation when the switching rates vary between 0.001 and 0.1 (Fig. S4) and between 0.001 and 0.5 (Fig. S5). Depending on the cooperation cost and selection pressure, there are still three main zones of the parameter space that correspond to the three evolutionary outcomes, three abundant/dominant genotypes: first-order cooperation, second-order cooperation, and pure defection. The exact borders between the zones move slightly, but all of our findings and the main message remain unchanged: second-order cooperators -non-cooperators that produce a large percentage of pure-cooperator offspring -thrive under a broad range of parameters for which pure cooperators are unable to invade.
A key difference between genotypic mutation and phenotypic switch is that the state of a phenotypic switch is often not inherited. So far we neglected this difference and all simulations have been conducted with genetic inheritance of the C/D alleles. To explore whether such inheritance makes a difference, we run a control simulation equivalent to the one presented on Figure 1 but without inheritance of the allele present at locus C/D: no matter the allele of the ancestor, the allele of a focal descendant will be C with probability M DC and D with probability 1 − M DC , where M DC is the mutation rate conferred by the allele H DC or L DC . Such model is much closer to the colicinogenic Escherichia coli and Salmonella typhimurium systems that we mentioned in the introduction. On Figure S3 we compare the outcome of this simulation (non-heritable cooperative behavior) with the simulation we performed on Figure 1 (heritable cooperative behavior). In what used to be the defection and second-order cooperation parameter zones when the switch was heritable, making the switch non-heritable does not change the outcome. This is expected because in this parameter zones, pure cooperators rarely reproduce but are mainly present as offsprings of second-order cooperators. However what was pure cooperation when the switch was heritable is now second-order cooperation, which makes sense because the cooperators can no longer form stable patches. These observations confirm that second-order cooperation can exist for both heritable and non-heritable types of cooperative behaviors.

Discussion
Our simulations have shown that the dominant evolved genotype in large spatially structured populations is often one that does not cooperate itself, but has a high probability of giving birth to offspring that do. We have termed the individuals bearing such genotypes the second-order cooperators. Using mathematical arguments and additional simulations, we show that high switching rate toward cooperation is a cooperative trait in and of itself, costly to the individual carrying it, but beneficial to others, thus being subject to the typical dilemma of cooperation. Second-order cooperation is maintained in a large portion of the parameter space, in spite of the pure cooperation behavior being extremely costly, or even effectively lethal, in the sense that pure cooperators almost never reproduce and mainly arise due to "production" by second-order cooperators. The maintenance is driven by the assortment between individuals bearing the allele that causes a high switching rate toward cooperation and the individuals who are cooperators. We can think of the second-order cooperation in two, not necessarily mutually exclusive ways, either as a division of labor or as a soma-germline distinction.
Both interpretations of second-order cooperation, division of labor, and soma-germline distinction, immediately point us to eusocial insects and the loss of reproductive behavior in worker casts, which has been extensively debated since mid-20th century (Hamilton 1964). While it is now well admitted that the spatial assortment between genetically similar individuals is central to the evolution of any cooperative behavior, including eusociality, the precise definition of this assortment is still largely debated, as well as the necessity and usefulness of explicitly representing it as a variable (usually called relatedness) in mathematical models (Wilson and Hölldobler 2005;Foster et al. 2006;Nowak et al. 2010;Abbot et al. 2011;Liao et al. 2015). A large part of the recent debate about relatedness is due to different authors using different definitions. For example while most modern authors define relatedness as a statistical measure of spatial assortment of similar individuals at a given locus -the one controlling the cooperative behavior - (Hamilton 1970), some do define it as a measure of similarity of several individuals across the whole genome (Wilson and Hölldobler 2005). In any case a high relatedness means that individuals interact more often with other individuals of the same genotype, for example that cooperators interact more often with other cooperators and thus benefit more from cooperation that non-cooperators. However, Fletcher and collaborators (2006) make the argument that relatedness is only one of many ways for cooperators to preferentially benefit from cooperation. Consequently, they insist that the assortment between the cooperative genotype and the effect of this genotype (e.g., worker phenotype) is broader than relatedness and is at the heart of the evolution of cooperation and eusociality. Correspondingly, the association that matters in our case is between the second-order cooperative allele (H DC ) and the effects of this allele, namely the presence of first-order cooperators, the individuals bearing allele C.
However, although we are insisting on allele and not wholegenome relatedness and are interpreting it as a spatial association, we should point out that this argument is not at all contradictory with inclusive fitness theory (Taylor and Frank 1996), despite the provocative statements made by Fletcher and collaborators. On the contrary, the association between our second-order cooperative allele and pure-cooperative phenotype falls into the extended inclusive fitness models developed by Queller (1985) and Grafen ( 1985). Additionally, it is similar to the spatial association between cooperators and public good in more classical settings: in our simulations pure cooperators are themselves a public good! Given that our populations are viscous, spatially structured, the effects of H DC allele, namely a higher presence of C individuals in neighborhood, are indeed directed toward other H DC bearing individuals more often than expected in a mixed population. Overall, although offering an important novel interpretation, our work is consistent and easily positioned within the historical and theoretical context of the research on the evolution of cooperation.
To confer upon the reader a more quantitative understanding of the effect of mutation rate and why it is creating a social dilemma, we suggest using a broader definition of fitness than the one classically used in game-theory inspired simulations. Instead of only considering the phenotype of a focal individual and thus its expected number of offspring, one should also consider the reproductive success of its descendants (Eshel 1973). Such definition of fitness makes sense because individual's genotype can affect the reproductive success of its the descendants independently from its own reproductive success, for example via mutation rate modifier alleles (Taddei et al. 1997). Put differently, when considering the long-term reproductive success of a lineage, the genotype of the ancestor contains more information than its phenotype alone, because it determines its position on the mutational landscape. This question of mutational landscape location as a selectable property of an individual is a much broader one, central to, for example, the theory of quasispecies (Eigen and Schuster 1977;Nowak 1992;Wilke et al. 2001;Colizzi and Hogeweg 2014).
Our work also directly connects with studies considering division of labor, another kind of social interactions between individuals. For example, Colizzi and Hogeweg (2014) present an RNA-world model where one "master sequence" encodes all the information necessary to produce (via mutations) a full ecosystem. Similarly, Rainey and Kerr (2010) suggest that non-cooperative types may be seen as a germline, whose ability to form a cooperative soma can be subject to selection. There is a common central theme between these two studies and the one we present here, which also corresponds to Darwin's "special difficulty": the ability to give birth to individuals with a certain phenotype can be a selectable property of a focal genotype that does not express the phenotype of interest.
Bacterial systems exhibit other types of division of labor, going beyond the soma-germline or queen-worker distinction, and relevant outside of the context of the evolution of cooperation. Some rely on phenotypic heterogeneity (Lopez et al. 2009), but some can rely on mutations (genotypic heterogeneity). It is the latter scenario that is observed in the case of the coexistence of several distinct lineages specialized in different resources, potentially including waste compounds produced by other lineages (Blount et al. 2008(Blount et al. , 2012Le Gac et al. 2012;Plucain 2014). We can identify two mechanisms that can maintain genotypic heterogeneity: the continuous generation of specialized genotypes by mutation from a master sequence (germline) and the stable coexistence of several specialized lineages. The question then becomes, which conditions may promote one versus the other mechanism. At first glance, the stable coexistence model may seem more plausible, at least under a reasonably low mutation rate, because it does not require second-order selection increasing the rate of a very specific mutation (contingency loci). However, the continuous generation model, exemplified by Colizzi and Hogeweg (2014), may be more relevant when mutation rates are very high, population bottlenecks are frequent, or selective sweeps are fast, making the coexistence of several lineages highly unlikely. All of the examples of stable coexistence we cited above indeed happened in Lenski's Long-Term Experimental Evolution (LTEE) lines, whose experimental setup is characterized by weak bottlenecks (1 in 100 each day, which corresponds to about 5 × 10 5 founder individuals). Interestingly, in the context of the diversification of a citrate-using subpopulation in LTEE experiments, Blount and collaborators (2008) also found that the switching rate toward cit+ genotype was increased while general mutation rate remained unchanged. This suggests that the two ways for genotypic heterogeneity maintenance are effectively not very distant or exclusive: while the two genotypes coexist for a long time, the ability of one to mutate into another may also be subject to selection.
Second-order cooperation can be supported not only by genetic mutations, which we focused on so far, but also by phenotypic switching. The primary differences between the two are (1) heritability, genetic mutations are inherited but phenotypic changes depend on likely weaker, epigenetic inheritance, and (2) occurrence, phenotypic switching is thought to happen at a higher rate than mutations. When the cooperative behavior is lethal, only the rate of occurrence may play a role. However, when cost of cooperation is lower and cooperators can coexist with non-cooperators, both may be important. For example, we can consider microbial species that occupy constantly changing natural habitats. We know that the majority of bacterial cooperative traits are only beneficial, and typically only expressed, in a fraction of the various environmental conditions bacteria experience (Kümmerli et al. 2009;Pai et al. 2012). A small fraction of the population may test for the potential benefits of cooperation by stochastically changing their behavior. While this could be done either via mutations or phenotypic switching, only mutations would create a stable, expanding subpopulation of cooperators. Phenotypic switching would result in less stable dynamics in the absence of a sensing mechanism, but may in turn be favored when the speed of change cannot be matched by mutations (Thattai and Van Oudenaarden 2004;Beaumont et al. 2009). Examples of both mechanisms exist in bacteria (Cascales et al. 2007;Ackermann et al. 2008;Hammerschmidt et al. 2014), but evaluating their relative importance requires further experimental and theoretical work.

Conclusions
In conclusion, we have shown that a non-standard type of cooperation, presenting some conceptual similarity with both division of labor and soma-germline differentiation, can evolve and be maintained even when cooperative behavior is very costly. This second-order cooperation, as we termed it, is characterized by dominant individuals that do not cooperate directly, but bear an allele that increases the chances of having pure cooperator offspring. We can interpret and better understand the evolutionary outcome of our simulations through spatial assortment. Indeed, we see high average spatial proximity (assortment) of secondorder cooperators to pure cooperators, making the former benefit from the presence of the latter, more so than the regular noncooperators who are typically located further away. Throughout the study we used synchronous, inelastic simulations based on prisoner's dilemma setup and maintain the generality of conclusions. Importantly, although our lattice simulations are not purely continuous, they have a spatial granularity as low as a single individual, and we do not need isolated subgroups to support our view. Our work highlights the importance of applying social evolution theory, for example kin selection, at a different level by considering not just the behavior of the individual but also of its offspring in the context of second-order selection. In terms of eusocial systems, rather than thinking about competition between workers and non-workers, we emphasize that it can be equally important to consider the competition between worker-generating and non-worker-generating genotypes. The evolutionary scenario we analyze is equally applicable to cases where the switching between cooperative and non-cooperative phenotypes happens via mutations or via regulation. Although some examples of secondorder selection in nature are known, including ones in the context of cooperation, we hope our results will motivate a closer experimental investigation of the prevalence and the range of conditions that can support second-order cooperation.

AUTHOR CONTRIBUTIONS
A.F. performed the simulations and analyzed the data. A.F. and D.M. wrote the manuscript. A.F., F.T., and D.M. designed the study, discussed, and interpreted the findings.

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's website: Figure S1: Allele frequencies for cooperation and mutation rate loci with randomized neighborhoods, when high mutation rate is 0.01. Figure S2: Invasion success of a rare H DC and L DC allele. Figure S3: Allele frequencies for phenotypic switch locus when high switching rate is 0.01. Figure S4: Allele frequencies for cooperation and mutation rate loci, when high mutation rate is 0.1. Figure S5: Allele frequencies for cooperation and mutation rate loci, when high mutation rate is 0.5. Figure S6: Allele frequencies for cooperation and mutation rate loci, when pure cooperators are sterile. Figure S7: Allele frequencies for cooperation and mutation rate loci, when simulations are asynchronous. Data S1