Why mussels stick together: spatial self-organization affects the evolution of cooperation

Cooperation with neighbours may be crucial for the persistence of populations in stressful environments. Yet, cooperation is often not evolutionarily stable, since non-cooperative individuals can reap the benefits of cooperation without having to pay the costs associated with cooperation. Here we show that active aggregation leading to self-organized spatial pattern formation can promote the evolution of cooperativeness. To this end, we study the effect of movement strategies on the evolution of cooperation in mussel beds. Mussels cooperate by attaching themselves to neighbours via byssal threads, thereby providing mutual protection. Using an individual-based model for mussel bed formation, we first demonstrate that the spatial pattern and the corresponding number of neighbours strongly depends on the movement strategies of the mussels. With an evolutionary model, we then show that this has important implications for the evolution of cooperation, since the evolved level of cooperativeness (the number of byssus threads produced) strongly depends on the number of neighbours and on the harshness and variability of the environment. Our results suggest that spatial aggregation, abundantly found in self-organized ecosystems, may promote the evolution of cooperation.

Abstract Cooperation with neighbours may be crucial for the persistence of populations in stressful environments. Yet, cooperation is often not evolutionarily stable, since noncooperative individuals can reap the benefits of cooperation without having to pay the costs associated with cooperation. Here we show that active aggregation leading to self-organized spatial pattern formation can promote the evolution of cooperativeness. To this end, we study the effect of movement strategies on the evolution of cooperation in mussel beds. Mussels cooperate by attaching themselves to neighbours via byssal threads, thereby providing mutual protection. Using an individual-based model for mussel bed formation, we first demonstrate that the spatial pattern and the corresponding number of neighbours strongly depends on the movement strategies of the mussels. With an evolutionary model, we then show that this has important implications for the evolution of cooperation, since the evolved level of cooperativeness (the number of byssus threads produced) strongly depends on the number of neighbours and on the harshness and variability of the environment. Our results suggest that spatial aggregation, abundantly found in self-organized ecosystems, may promote the evolution of cooperation.

Introduction
Cooperation between neighbouring individuals is often essential for survival in stressful environments (Bertness and Callaway 1994;Callaway and Walker 1997;Holmgren et al. 1997;Stachowicz 2001). Organisms ameliorate their environment locally, for instance by providing shade or by drawing moisture and nutrients towards themselves and close neighbours (Schlesinger et al. 1996;Aguiar and Sala 1999), which allows others to survive in an otherwise hostile world. To what extent cooperation evolves in a population depends on the nature and intensity of interactions between individuals (Hamilton 1963;Axelrod and Hamilton 1981;Nowak and May 1992;Rainey and Rainey 2003;Doebeli and Hauert 2005;Santos and Pacheco 2005;Ohtsuki et al. 2006;Santos et al. 2006;West et al. 2007;Masuda 2007;Van Dyken and Wade 2012). The number of cooperating individuals an organism interacts with likely determines the effectiveness of its cooperation strategy and may affect the degree of cooperativeness that evolves within a population (Vainstein and Arenzon 2001;Zhang et al. 2005;Ohtsuki et al. 2006;Hui and McGeoch 2007).
Even if cooperation is profitable for all interacting individuals, it is intrinsically unstable when recipients can reap the benefits of cooperation without helping others in return. Such a social dilemma can be solved to a certain extent by spatial population structure (Axelrod and Hamilton 1981;Nowak and May 1992;West et al. 2002). Especially in highly viscous populations, where cooperative traits are transferred locally, the effects of network structure and neighborhood size have been studied extensively (Pfeiffer and Bonhoeffer 2003;Rainey and Rainey 2003;Santos and Pacheco 2005;Ohtsuki et al. 2006;Santos et al. 2006;Masuda 2007). Yet, little is known about the evolution of local cooperation in species that disperse their offspring over a wide range, but interact locally. This is, in particular, the case for widely dispersing organisms that upon settlement, move into a self-organized spatial structure.
Systems as diverse as mussel beds, coral reefs, marsh tussocks, tidal wetlands, peat lands, arid ecosystems, and ribbon forests are highly structured in space due to the interplay of local facilitation and long-range inhibition (Klausmeier 1999;Mistr and Bercovici 2003;Rietkerk et al. 2004a, b;van de Koppel et al. 2005van de Koppel et al. , 2008van de Koppel and Crain 2006;Rietkerk and Van de Koppel 2008;Eppinga et al. 2009). In these systems, the number of potentially cooperating neighbours depends on the spatial scale and distribution pattern of the population. In many systems, the spatial pattern results from the active movement of organisms (Theraulaz et al. 2003;Jeanson et al. 2005;van de Koppel et al. 2008;de Jager et al. 2011;Hemelrijk and Hildenbrandt 2012). Accordingly, the movement strategies of these organisms can indirectly affect the number of neighbours an individual will encounter. In situations where costs and benefits of facilitation depend on the availability and density of local neighbours, the movement strategy therefore affects the evolution of facilitation.
An example of active pattern formation can be found in intertidal mussel beds. Mussels self-organize into large-scale labyrinth-like patterns (van de Koppel et al. 2005. They use their foot to aggregate into a group of conspecifics after wide dispersion by the currents during the larval stage (Geesteranus 1942). Because mussels are well-mixed during their larval stage, relatedness between neighbouring individuals is, on average, equal to the relatedness between distant individuals within the same mussel bed (Ferguson et al. 2013). When aggregated, mussels facilitate each other by attaching byssus threads (a glue-like substance made of protein strands, which are costly to produce; Eckroat and Steele 1993) to the shells of conspecifics that are within reach. These attachments decrease the risks of dislodgement and predation for both the attaching mussel and the one receiving the byssus thread Scheibling 2001, 2002). Mussels that are sufficiently affixed by neighbours do not need to create attachments themselves and can therefore avoid the costs of producing byssus threads. Through active aggregation into mussel clumps with various densities, mussels can modify the number of neighbours within their attachment range. By self-organizing into the labyrinth-like patterns that are characteristic for intertidal mussel beds, mussels attain an intermediate number of neighbours, which lies between the few neighbours that are within attachment distance in scattered distributions and many neighbours in dense mussel clumps.
In this paper, we investigate how spatial patterns affect the evolution of cooperativeness in self-organized mussel beds. For this purpose, three questions regarding cooperation in mussel beds will be addressed. First, we investigate how the aggregation strategy of mussels affects the spatial pattern and, in particular, the number of neighbours available for cooperation. Aggregation in mussels typically leads to the formation of a spatial pattern consisting of regularly spaced strings and clumps (van de Koppel et al. 2005. With an individual-based model (IBM), we investigate how the number of neighbours a mussel experiences is related to this self-organized pattern. Second, we examine, by means of an adaptive dynamics approach, how the number of neighbours affects the tendency to attach costly byssus threads to neighbours, which we interpret as the cooperativeness of an individual. Building on the fundamental assumption that the spatial pattern relates to the average number of neighbours that a mussel can attach its byssus threads to, investigating how the number of neighbours affects the evolution of the attachment tendency of mussels gives us insight into whether and how aggregation strategies promote or hamper cooperation. Third, we study the effect of harshness of the environment. It is likely that this affects the evolution of cooperation, since for mussels, survival under stressful conditions depends on how well they are attached to their neighbours. Furthermore, we take into account that environmental stress likely differs substantially between generations, which may further affect the evolution of cooperativeness.

An individual-based model of self-organized patterning
We modelled the effect of individual aggregation strategies on the formation of mussel beds with an IBM. As the self-organized pattern in mussel beds is a compromise between reducing wave stress and predation risk (requiring dense aggregations) on the one hand and minimizing food competition (requiring low densities on a larger spatial scale) on the other (van de Koppel et al. 2005, mussels move around until they find a location where local mussel densities are sufficiently high and densities at a larger scale are low enough to permanently establish in the mussel bed. We developed an IBM that describes pattern formation in mussels by relating the chance of movement to the density of the mussels at two spatial scales, i.e. within a short-distance of 2 cm and a long distance of 7.5 cm, following de Jager et al. (2011). We consider 1600 circular individuals with a diameter of 1 cm that are initially spread homogeneously on a 30 9 30 cm surface. In each of the 500 time steps within a simulation, all individuals get a chance to move in random order. Whether a mussel moves or not depends on the density of mussels within the local attachment range of 1.1 cm ø (i.e. the 'local density') and the density of mussels within the larger, 3.3 cm ø competition range (i.e. the 'long-range density'); a mussel moves when the local density is lower than a certain settlement threshold (which we will vary below) and/or when the long-range density is higher than 0.7 individuals/cm 2 . These parameter values were estimated using a regression analysis of experimental data de Jager et al. 2011). We modelled the movement of individuals in correspondance to natural mussel movements, using a heavy-tailed step length distribution (a Lévy walk with l = 2; de Jager et al. 2011), where steps are made in random directions and their lengths are drawn from a power law distribution. A mussel ends its step prematurely when it encounters a conspecific (de Jager et al. 2014). In our model, mussels cooperate after (and not during) pattern formation; therefore the attachment of byssus threads does not impair mussel movement. To examine the relation between the number of neighbours within the facilitation range and the spatial structure that emerges in the self-organized mussel bed, we simulated mussel bed formation for a range of settlement thresholds, e.g. the minimum mussel density required for local aggregation. We plotted the emergent spatial patterns and calculated the average number of neighbours ±SE within attachment range for each simulation.

A model for the evolution of between-mussel cooperation
To investigate the evolution of cooperation, we make two plausible assumptions on how the survival probability and the fecundity of a mussel are affected by its attachement tendency A and on the number n of neighbours within attachment distance. The attachment tendency A (0 B A B n) corresponds to the average number of byssus threads produced by a mussel and attached to its neighbours. Mussels, however, do not only make attachments themselves, but also receive attachments from other mussels. Hence, the total number of attached neighbours N depends on both a mussel's own production of byssus threads (A) and on the number of attachments produced by its neighbours. A mussel can be attached to a neighbour by its own byssus thread, by the byssal attachment of its neighbour, or by both; it stays disconnected from the neighbour if both do not attach to one another. Thus, we can calculate the probability that two mussels are attached as 1 minus the probability that they remain disconnected. Given that a mussel has n neighbours, an attachment tendency A, and neighbours with an attachment tendency A 0 , the expected total number of attached neighbours is given by We consider this total number of attached neigbours to be an important determinant of an individual's survival probability.
We assume a nonlinear relation between the number of attached neighbours and individual survival probability (Archetti and Scheuring 2011), that survival is high when a mussel is attached to many neighbours and is much lower when a mussel has only few attached neighbours: Here E is the number of attached neighbours needed for the survival chance to be 50% and k determines the steepness of the logistic, S-shaped function (Fig. 1). Throughout, we will assume that survival for mussels attached to zero neighbours is 1% (S 0 (0) = 0.01). This imposes a constraint on the parameters k and E, essentially reducing the number of parameters to one. We further assume that the production and attachment of byssus threads has fecundity costs and consider a linear relation between fecundity and the average number of byssus threads produced: Here c denotes the costs per cooperation with a neighbour (Nicastro et al. 2009). Note that the production of a byssus thread can either be directly beneficial as well as costly, depending on the number of neighbours already attached (an additional byssus thread will only slightly increase survival when N(A, A 0 ) is large, and costs will likely outweigh benefits) and the encountered level of environmental stress (E). For simplicity, we do not consider responsive or conditional strategies in our model.
To study the evolution of the attachment tendency, we use an adaptive dynamics approach (Geritz et al. 1998). To this end, consider a monomorphic resident population with attachment tendency A 0 , in which a mutant with strategy A arises. Whether this mutant invades the resident population depends on its relative fitness (W). For simplicity, we assume that the individuals in our model are semelparous. Then fitness corresponds to the product of the probability to survive (S) until reproduction and expected fecundity (F). Hence, the relative fitness of a mutant with attachment tendency A is given by: If W(A, A 0 ) [ 1, the mutant genotype has larger fitness than the resident genotype and can increase in relative frequency. Assuming asexual reproduction and mutations of small effect, the invasion of a mutant when rare typically guarantees that the mutant will spread to fixation, hence replacing the former resident (Geritz et al. 1998). Through a series of consecutive gene-substitution events, the attachment tendency will evolve to an evolutionarily singular strategy A* (Dercole and Rinaldi 2008). Such a strategy is evolutionarily stable if no mutant strategy can invade a population of individuals using strategy A*. An Evolutionarily Singular Strategy A* is convergence stable if those mutants successfully invade a given resident strategy A 0 that is closer to A* (Geritz et al. 1998).
The parameter E in Eq. 3 represents environmental conditions, such as wave stress and predation risk. In harsh environments, E will take on a larger value than in benign environments. We will examine the evolution of attachment for a range of environmental conditions. Furthermore, environmental conditions are likely to vary between generations. Hence, we will also investigate the effect of alternating environments on the evolution of cooperation.

Spatial patterning relates to number of neighbours
As a first step, we demonstrate that the aggregation strategy of mussels strongly affects their spatial distribution as well as the number of neighbours a mussel can interact with. To this end, we systematically changed the settlement threshold of the mussels in a population. Simulations of our individual-based model reveal that a scattered distribution results when the settlement threshold is low, that a labyrinth-like pattern emerges when the settlement threshold is intermediate, and that dense clumps are formed when the settlement threshold is high (Fig. 2 top). The average number of neighbours increases with the degree of aggregation (Fig. 2 bottom). Using different average mussel densities in the simulations results in a similar pattern, though the number of neighbours increases with the overall density (Fig. 3). For the remainder of this paper, we will use the evolutionary model described above, thus making the assumption that the number of neighbours (n) represents a certain degree of aggregation (for a given mussel density). We will use the following number of neighbours (n) to represent the different spatial structures: n = 6 for scattered distributions, n = 8 for labyrinth-like patterns, and n = 12 for dense mussel clumps. Because natural mussel beds are often labyrinth-like, we are specifically interested in how an intermediate number of neighbours (n = 8) affects the evolution of the attachment tendency A.

Evolution of the attachment tendency A
For three different environmental conditions [benign (E = 2), moderate (E = 6), and stressful (E = 10)], Figure 4a shows how the evolutionarily stable attachment strategy A* depends on the number of neighbours within attachment range n. Interestingly, the number of neighbours for which investment in cooperation is maximized increases with the level of environmental stress. Investment in attachment peaks at different numbers of neighbours for different levels of environmental stress. In Fig. 4a, A* first increases more or less linearly with n before levelling off. When the number of neighbours within the attachment range is low, individuals attach themselves to virtually all their neighbours (A* & n for small values of n). With increasing n, the mussels attach themselves to a smaller and smaller percentage of neighbours; in fact the number of byssus produced declines with n when many neighbours are available. Depending on environmental conditions, cooperativeness (=the number A* of byssus threads produced) is maximal at a value of n that corresponds to a labyrinth-like pattern or a dense mussel clump.

Changing environmental stress levels
Because mussels disperse over a wide range as larvae before settling on a mussel bed, environmental conditions are most likely different between generations. Adaptation of between-mussel cooperation to a particular stress level is therefore difficult and evolution of cooperation becomes more challenging than described above. We investigate the Fig. 3 The average number of neighbours increases with overall mussel density. For each mussel density, mussels in labyrinth-like patterns have, on average, more neighbours than those in scattered distributions, but less than those in dense mussel clumps. For low densities, mussels were unable to aggregate into labyrinths or clumps. Average mussel density was calculated as the total number of simulated individuals (M) per surface area (M/0.30 2 m) Evol Ecol (2017) 31:547-558 553 robustness of our results above to variability in the environmental conditions that are encountered during the adaptation process. In Fig. 4b, we considered the three situations where the environmental stress level a generation encounters is variable; drawn from a random distribution (l = 6) with low (r = 1), intermediate (r = 3), and high (r = 5) variation in stress. When variation in E is high, the evolutionarily stable attachment tendency is very low for all n. With a mean stress level l = 6, only at low variation in environmental stress do we find a hump-shaped relation between the number of neighbours and the average number of attachments a mussel produces. This confirms the results we obtained in the absence of environmental variation between generations (Fig. 4a). Fig. 4 a The relation between the number of neighbours and investment in the number of attachments created to neighbouring individuals is hump-shaped and varies with the level of environmental stress. Note that the number of attachments created can never be larger than n. b Evolution of attachment tendency when environmental conditions differ between generations and vary according to a normal distribution. The solid line indicates the case where environmental stress is normally distributed with little variance (l = 6, r = 1); variance is increased for the two dashed lines (r = 3 and r = 5, respectively). The grey areas indicate the error margins (average ± SE) of 10 simulation runs with the same parameter settings Fig. 5 The average level of environmental stress (l) and the inter-generational variation in stress (r) determine which spatial population structure (scattered, labyrinth, or clumped) maximizes the investment in between-mussel cooperation. The figures were generated by interpolating the results of 21 9 21 simulation runs with different combinations of average environmental stress and inter-generational variation in stress. Simulations were run with 1000 (a), 1600 (b), and 2000 (c) individuals The effect of inter-generational variation in stress becomes more apparent in Fig. 5, where we show the relation between mean environmental stress, inter-generational variation in stress, and the spatial structure that results in the highest number of attachments produced. Simulations with different overall mussel densities show similar patterns: with increasing average stress and/or variance in stress, the spatial structure that maximizes investment into cooperativeness shifts from a scattered distribution, to labyrinth-like patterns, to dense mussel clumps (Fig. 5). These results suggest that intermediate stress levels together with low to intermediate inter-generational variation in stress can cause mussels in labyrinth-like patterns to evolve higher investment in byssal attachments than mussels in scattered distributions and dense mussel clumps. In other words, the labyrinth-like pattern that we observe in nature only promotes between-mussel cooperation better than other spatial structures under a limited set of environmental conditions.

Discussion
Cooperation is often a necessity for survival in harsh environments and is therefore found in many species. Organisms utilize a multitude of supporting traits and behaviours, such as local dispersal, reciprocity, and punishment, to maintain high levels of cooperation (West et al. 2007). Our theoretical analysis reveals that in intertidal mussels, movement into spatial aggregations stimulates the evolution of cooperation. Because mussels benefit from any attachment of byssus threads with neighbouring individuals, some degree of betweenmussel cooperation evolves in any type of mussel bed, irrespective of the number of neighbours. However, our analysis shows that the number of neighbours that maximizes investment in cooperation depends on environmental conditions and overall mussel density. In low stress environments with little inter-generational variation in stress, aggregating in scattered distributions maximizes investment in cooperation. In contrast, investment in cooperation is maximized when mussels aggregate in dense clumps in high stress environments with considerable variation in stress between generations. Yet, aggregating in labyrinth-like patterns, which mussels do in natural mussel beds, only maximizes investment in byssal attachments in a small range of environments with intermediate stress levels and corresponding inter-generational variation in stress. Based on our results and those of others (Ohtsuki et al. 2006;Santos et al. 2006;Masuda 2007), we can conclude that forming spatial aggregations can substantially influence the degree of cooperativeness that evolves in a population.
For simplicity, we did not take the correlation between environmental stress and food availability into account. In most intertidal ecosystems, an extensive range of environmental conditions can be encountered at any time, from very benign habitats that also provide little food, to very hash conditions where food is often abundant. Mussel offspring is likely to reach all of these habitats, as is wittnessed by the high availability of mussel spat on artificial settlement structures. This implies that the offspring of any mussels can spread itself over different habitats where a more harsh environment implies a better food supply. Further research may show whether the inclusion of this relationship between environmental stress and food availability will give different results. It is likely that the levels of cooperation that are found in real-world mussels reflects an adaptation to the habitat where they can generate the highest number of offspring, taking into account the availability of the habitat in the overall area.
We furthermore adopted a number of simplifying assumptions that do not precisely reflect the conditions that mussels, or any real-world organism, would encounter. In mussels, reproductive output per unit of biomass increases with age, as growth takes an ever smaller part of energy. Under most circumstances, our simplification of semelparity has little consequences, yet it might become important in temporally variable environments. We assumed a fixed self-organizing behavior within each and throughout generations; in each simulation of our IBM, all individuals used the same set of rules, including the settlement threshold, to move into a spatial pattern. This is an unrealistic assumption for several reasons. For example, generations are likely to differ in initial overall density; a scattered population in a dense mussel bed will result in a higher number of neighbours within attachment distance than in less dense but patterned beds. In this case, mussels in our model will never stop moving and hence never attach to any neighbours, because the long-range mussel density remains too high. Furthermore, individuals might differ in their self-organizing strategy; though some are aggregating in dense clumps, others may be strategically moving away from dense mussel clusters. A further simplification is that we only examined one aspect of spatial patterning on mussel survival: the effect of the number of direct neighbours. Irrespective of the overall mussel density, the number of neighbours is lowest for a scattered distribution, intermediate for a labyrinth-like pattern, and highest in dense clumps. Still, depending on the overall mussel density, an equal number of neighbours, and hence a similar degree of cooperation, can be achieved in all three patterns. It is conceivable that not only the number of neighbours, but also the spatial patterning of the neighbourhood is of importance for the evolution of cooperation. Spatial population structure plays an important role in mussel beds, as it defines not only the number of primary connections, but also secondary and tertiary connections between mussels, which bond many mussels into a single clump. Production of clumps generates an additional selection pressure, as larger clumps are less likely to become dislodged by wave stress. Though we do not consider the effect of spatial patterning on group size and higherorder selection processes in the current paper, we do analyse these effects in a separate study (de Jager 2015).
Our study demonstrates that active self-organization can have substantial consequences for the degree of cooperation between neighbours that evolves in a population. Inversely, self-organized spatial patterns have been described in a wide range of ecosystems, and many of these studies highlight the importance of cooperative interactions for the formation of these spatial patterns. In patterned arid bushlands, for instance, plants promote the infitration of water into the soil, facilitating other plants (Klausmeier 1999). This highlights the potential importance of feedback interaction between pattern formation processes on the one hand, and cooperation on the other. Whereas studies on the evolution of multicellularity have shown that feedback between pattern formation and evolution of cooperation/division of labour can drive evolution of unicellular to multicellular organisms (Pfeiffer and Bonhoeffer 2003;Rainey and Rainey 2003;Ratcliff et al. 2012), our study system may provide a suitable template to investigate such feedback in populations of nonrelated individuals. So far, the evolution of cooperative interactions other than aggregation and the pattern forming characteristics of organisms, such as their aggregative behavior, have been studied in isolation. Although our conclusions-that evolution of cooperation depends on spatial aggregation within the population-can be drawn without the explicit inclusion of joint evolution, the joint evolution of pattern forming and cooperative traits is a promising subject for further investigation.