Speciation, Diversification, and Coexistence of Sessile Species That Compete for Space

Speciation, diversification, and competition between species challenge the stability of complex ecosystems. Laboratory experiments often focus on one or two species competing under conditions where they may grow exponentially. Field studies, in contrast, emphasize multi-species communities characterized by many types of ecological interactions. A general problem is to understand conditions that support a dynamically maintained coexistence of many species in an ecosystem over a long time span. In the present paper we propose a lattice model of multiple competing and evolving sessile species. When allowing the interspecies interactions to mutate, we obtain coexistence of many species in a complex ecosystem, provided that there is a cost for each interaction. The diversity reached by the model incorporating speciation is found to be substantially higher than in the case when entirely new species appear due to immigration from outside of the considered ecosystem. The species self-organize their spatial distribution through competitive interactions to create many patches, implicitly protecting each other from competitively superior species, and speciation in each patch leads the system to high diversity. We also show that species that exist a long time tend to have a relatively small population, as this allows them to avoid encounter with competitive invaders.


Introduction
Biological organisms cooperate and compete with each other forming complicated ecological systems with an intriguing ability to sustain themselves over long time-spans. In fact, this stability is not easy to understand as the interplay between exponential growth and competition produces an inherently unstable state. Accordingly, ecosystems consisting of more than a few species should tend to collapse into a low diversity state. In a seminal paper R. May pinpointed that this instability is weakened by reducing the number of interactions in the ecosystem [1]. One way to reduce interactions as well as exponential growth and competition is to introduce spatial segregation [2][3][4][5][6][7][8][9].
Many ecosystems consist of multiple interacting species that may form niches for each other [10][11][12], exemplified by the concept of keystone species [13,[13][14][15]]. As a model for selforganized niche formation, a simplified description of competing lichen species on the two-dimensional surface of a rock has been introduced in Refs. [16,17]. This model considers the spreading and competition of mutually exclusive species on a two-dimensional lattice. The competitive interactions between species are assigned randomly. Each lattice site contains a maximum of one species which attempts to colonize the neighboring sites. The invasion is possible only if the ecological interaction between the invader and the species that occupies the neighboring site allows it. In addition, new species are occasionally introduced at random positions, leading to a state of high diversity, provided that the likelihood for interactions is low [16,17]. If the model is simulated in the well mixed situation by allowing interactions between spatially separated species, the high-diversity state collapses [16]; the space is essential for the high-diversity state. The short range interactions are one of the important sources for the high-diversity [7]. Another important factor is transient cyclic invasion that generates patches of isolated niches when it collapses [17].
In the present work, we investigate a two-dimensional evolution model for sessile species, where new species are not introduced from outside, but instead evolve from the already existing species. The major result observed for the model investigated in Refs. [16,17] is found also now, i.e., the high diversity state is stable when the species invasive interactions are sufficiently inhibited. The model allows us to combine allopatric and sympatric speciation through a self-organized segregation of species into isolated patches. The mutations in segregated patches allow neutral evolution and lead to evolutionary divergence of the properties of the spatially separated species.

Model
We investigate a model representing the evolution of an ecosystem that initially consists of D~1 species only, occupying one randomly chosen lattice site of a L|L square lattice with periodic boundary conditions. For simplicity, we assume that a species is characterized by its ecological interactions, expressed by the interaction matrix C as explained below. Namely, the only phenotype we focus on, is the interspecies interaction; we ignore the difference between the genotype and phenotype. The differences in the ecological interactions do not necessary indicate the difference in, for example, interbreeding ability. Therefore the word strains instead of species may be sometimes more suitable when the difference between diverging lineages is small, though after long time the accumulation of divergent mutations tend to promote the origin of widely different species. In the following we use the term species only for simplicity.
At each update a site i is selected randomly among the occupied sites. With probability a N the species s(i) mutates to change its ecological interactions. If no mutation takes place, a site j is selected randomly among the four nearest neighbors of site i. If the species s(i) can invade the species s(j), C½s(i),s(j)~1, then the corresponding invasion takes place, i.e., site j is updated by setting s(j)~s(i). Here C is a matrix that represents possible interactions between the species, which remain fixed once they are introduced. It is assumed that if s(i)~s(j) then C½s(i),s(j)~0 and that an empty site can always be invaded, C½s(i),0~1. One Monte Carlo time step is defined as L 2 repetitions of the described updates.
In the case of mutation the interaction matrix elements of the new species s Ã coming from the ancestor s(i) are assigned according to the following rules: 1. for any existing species s' one initializes C½s Ã ,s'~C½s(i),s' and C½s',s Ã ~C½s',s(i); 2. one randomly selected element of C½s Ã ,s' or C½s',s Ã is set to 1, i.e., if the value was 0 then now the mutant s Ã can invade a species s' or can be invaded by a species s' which had no interaction link with the ancestor s(i), respectively; if the value was already 1 then nothing changes; 3. for each s' the elements C½s Ã ,s'~1 and C½s',s Ã ~1 are set to 0 with probability E, representing a cost in maintaining interactions; 4. finally, we set C½s Ã ,s(i)~1 and C½s(i),s Ã ~0, i.e., it is assumed that the mutant s Ã can invade its ancestor s(i) but the ancestor cannot invade its mutant.
In short, the new species s Ã inherits most of the phenotypical features of its ascendant s(i), represented in the competition network, with small random modifications. This is in contrast to the previously studied model in Refs. [16,17], where the interactions for a new invading species was assigned randomly with a pre-determined interaction probability c.
The rule 4 that new species can always invade its ancestor while the ancestor cannot invade its mutant is motivated by the following: One could consider that one site is already a collection of many individuals. Then, C½s Ã ,s(i)~1, C½s(i),s Ã ~0 is the combination that the new species will succeed to fixate. If C½s Ã ,s(i)~0, s Ã is unlikely to take over already existing s(i) in the site. Even if C½s Ã ,s(i)~1, if C½s(i),s Ã ~1 also holds then it is expected to be hard for the species s Ã to win over majority species s(i). In addition, if we open for the stand-off relations, namely C½s Ã ,s(i)~0 and C½s(i),s Ã ~0, then this will create a lot of small static dusts of new species in the ancestor species region. Since the separated patches are crucial for higher species diversity as we have already seen in the original model [16,17], having this additional possibility is expected to increase the diversity further more, but such an effect may be too artificial.
The time scale of mutation is expected to be much slower than the time scale of interspecies interactions in an ecosystem. Therefore, in some of our simulations we investigate the quasistatic limit, which effectively is equivalent to the limit a N ?0. This approximation focuses on the evolution of a new species, whereas its effect on population redistribution is speeding up cyclic interactions to re-establish a representative frozen state [16] before the next mutation. New species appear through mutations only after all the activity of invasions has died out, because when multiple species compete for a finite region, the stochasticity of the dynamics eventually makes one of the species to take over. In most of the cases this process takes quite a short time, about L steps, which is the time scale for a front to sweep the whole system. However, there are occasional events where several species compete dynamically for the same area for a long time period, often due to a short cyclic relationship. To save the computation time, we shorten the long-lasting competition by temporally preventing one of the randomly chosen active species to invade any other species after typically 100|L time steps. It should be noted that the one species state is an absorbing state in this limit; therefore we use the high-diversity state obtained by a finite a N simulation as the initial condition for the quasi-static simulation [16,17]. In order to reduce the computation time, the simulations with small values of a N , including the quasi-static simulations, were performed by the event-driven type algorithm, where the possible events are listed and time to the next event was drawn accordingly from the exponential distribution. We have verified that this method gives statistically the same results as the described random sequential updates. Typically the model is simulated for a long time before a reliable analysis of steady state properties can be made. The subsequent section analyzes aspects of this steady state dynamics.

Time evolution of the system
In order to understand the dynamics of the system, let us start by investigating the time evolution of a stochastic realization of the system. Figure 1A presents snapshots of the system at different times after the steady state has been reached. Figures 1B and 1C present the time evolution of population size N of the species marked in red in Fig. 1A and of its offspring species during a shorter and a longer time interval, respectively. The rest of the species in the snapshots are shown with various shade of blues, except for one of the off-springs (yellow) and the species that invade the yellow species (green).
By comparing the snapshots 1 and 2 in Fig. 1A, we can observe the expansion of the species marked in red in the territory that was initially occupied by one of the species marked in blue. This expansion corresponds to a step-like change in its population size, N, in Fig. 1B. Subsequently, a long silent co-existence of the two species takes place (compare snapshots 2 and 3 in Fig. 1A), where the population size, N, of the red species remains constant (Fig. 1B). The mutation rate a N~5 |10 {6 corresponds to one mutation per 500 time-steps among about 500 lattice sites of the red population. The change in population size coincides with mutations in red species with this time scale; there occurs a mutation of the red species with a speciation into a yellow descendant (see snapshot 3 in Fig. 1A). The yellow species does not expand into the territory of the blue species, but it invades its ancestor region (compare snapshots 3 and 4 in Fig. 1A). However, the expansion of the green species separates the red species and its yellow offspring (see snapshot 5 in Fig. 1A) at time t&4900. Finally, the green species brings the yellow one into extinction, letting at the same time the red species to survive (see snapshot 6 in Fig. 1A). Figure 1C illustrates the fate of the red species on a longer timescale. The dynamics is persistently punctuated [18], with small time intervals of rapid changes followed by long periods where the population size, N, is constant. In terms of the fitness concept, the growth rate l of the population, is associated to the fitness F as F~log (l), and F deviates from 0 in an erratic way only for short time intervals. The inserts below the population size, N, in Figs. 1C and 1B show the growth rate, l, of the population, which can be negative when the species is exposed to its competitively superior species. Whereas the measurement of l in a short time scale naturally reflects the contemporary environment, the assignment of the fitness in a long time scale is meaningless. Instead, the survival of a species in a longer time scale depends on the extent to which the species hedges its option by occupying separate patches. Figure 2 examines the larger scale dynamics of the system with size L~100 for a value of E leading to the situation when the system may persist at both low and high diversity. Figure 2A shows the time dependence of the diversity, D, and the number of patches, P, in the system after the steady state has been reached. One can observe the switching between a state with low diversity and highly varying patchiness, and a state with high diversity and high patchiness. This may resemble the bi-stability found earlier in the model of Refs. [16,17] for the interaction density c&0:07 (c is the probability that an element in C-matrix takes the value 1), but in the present case c changes with time (see below). There can also be a finite size effect, regarding that the previous model required Lw150 to have a stable high diversity state [16,17].
At any time, the system may be characterized by the number of species, D, and by the interaction density, c, that can be calculated from the contemporary interaction matrix C at that time point as follows: c~P D s1~1 P D s2~1 C(s 1 ,s 2 )=½N(N{1) (remember that the diagonal element of the interaction matrix, C, is zero). Figure 2B examines c, which takes the value c&0:03 in the high diversity state and is oscillating between 0 and 1 in the low diversity state. Figure 2B therefore illustrates a positive feedback between the high diversity and the difficulty in maintaining potential interactions with many species. For the high diversity state, the value of c can be estimated quite accurately (data not shown) from the steady state condition by adding and removing links assuming that the matrix is represented by c: Steady state properties The steady state properties of the system are presented in Fig. 3 for a lattice size L~200; all panels in Fig. 3 examine three different values of the mutation rate a N , in order to investigate the limit a N ?0. To see the transition between the low and high diversity states, we study in Figs. 3A and 3B, respectively, the average diversity, SDT, and the average interaction density, ScT,  versus E (the averages are taken over a long time after the steady state has been reached). As one can see from Fig. 3A, in the quasistatic limit, a N ?0, there is a sharp transition from low to high diversity state between the values E~0:01 and E~0:015. The figure indicates also that higher speciation rates, a N , are associated to the high diversity state that can be maintained at substantially lower interaction penalties E. It was confirmed that the well-mixed version of the model did not give a stable high-diversity state, namely the space is necessary (data not shown).
The comparison of Figs. 3A and 3B shows that all high diversity states are associated to c values below a critical value that is about 0:06, apparent also from the function SDT versus ScT in Fig. 3C.  For comparison, Fig. 3C shows also the behavior of SDT versus c in the non-evolutionary model of Ref. [16], recapitulating the critical c value for this version of the model. In general, the qualitative first order transition between the low and high diversity states is found to be robust against the rule that species interaction is ''inherited'', i.e., copied to new species with only small changes. The main quantitative difference with the non-evolutionary model of Ref. [16] is that the present version of the model leads to approximately three times larger value of D in the high diversity state (see Fig. 3C).
In order to gain more understanding about the system, it is useful to investigate also the Hamming distance between the species m and n, characterizing their similarity, The Hamming distance in the high diversity state is very close to the expected value for randomly assigned C-matrix with the interaction density ScT, when averaged over the neighboring species (H neighbor ) and over all pairs of species in the system (H all pair ), as shown in Table 1. This result indicates that species mutate many times neutrally before spreading substantially in space. This local accumulation of mutations suppresses spatial correlations between the phenotypes of neighboring species.

Following lineages
The proposed evolution model allows one to follow descendants of individual species and therefore to examine to what extent the change in population size is correlated with the diversity and the survival of a lineage. Figure 4 shows one case history, following a fairly successful species that for a limited time period has more than ten different offspring species occupying up to 4000 lattice sites on a 200|200 lattice. Whereas the number of offspring   The same data as in panels A and B plotted as SDT versus the measured interaction density ScT: Note that ScT decreases to the right. In the case of a N~1 0 {7 and in the quasi-static simulations we started from the high diversity state obtained for a N~1 0 {6 and E~0:04 and simulated the relaxation to a constant value. For comparison, the results for the model without heredity, studied in Ref. [16], are also plotted in the quasi-static limit. For a smaller system size, e.g., L~100, the high diversity states shifts toward E?1 as a N ?0, indicating that the high diversity state is not stable in the quasi-static limit for a too small system (data not shown). doi:10.1371/journal.pone.0096665.g003 species develops quite gradually, one can notice that the population changes are intermittent. Remarkably, the lineage mostly persists in a stage with relatively few individuals. Figure 5 systematizes the population statistics of species and their families as a function of their life-time (a N~1 0 {7 , E~0:02, and L~200). Interestingly, the species that exist only for a rather short time (10 3 -10 4 time units) are those whose population sizes are the largest. In contrast, the species that exist for sufficiently long time to allow speciation (blue circles) mostly have small population sizes. We emphasize that also the maximum population of a species is declining when its total existence time becomes long (data not shown). Instead, the red circles depict the average population size, N, of the whole lineage starting with a single species. Again, one can observe a decline in population as survival time is longer, a pattern that is also consistent with the case story plotted in Fig. 4.

Discussion
In the present work we investigated a model of competition between sessile species that evolve. As we have defined a species in terms of interaction options then the evolution is represented by small changes (mutations) in these links. The focus of the proposed model is speciation, where a population of a species is spatially disconnected by interspecific competition, and this disjunction facilitates subsequent divergence by allowing each patch to mutate and develop differently. This speciation scenario resembles to some extent the allopatric speciation for evolution of plants [19,20]. However, here it is determined by other species, and not by mountain rages or other non-biological factors. As far as we know, this is a new proposal of a mechanism to create long-term isolated populations (meta-populations). The closest ecological data related to spatially driven speciation is the studies on how the dynamical change of geological environment affects the speciation events [21][22][23]. These studies consider the situation where the geological changes are due to external causes such as climate change, but it will be interesting to include the possibility that the change of species spatial distribution is caused by the interspecies interaction.
We have observed that the population in each patch evolves in a series of replacements where new mutants replace the old ones. Most of these replacements are phenotypically neutral in the given environment of the patch, allowing large scale evolutionary meandering in analogy to neutral evolution [24][25][26]. As a consequence of this extensive neutral development within each path, the species in neighbor patches are typically as different from each other as any random pair of species. When the inhabitant of a patch finally becomes capable to invade a neighbor, it emerges as a species that is entirely different from the species that originally established the patch. Therefore the model behaves to a large extent as the immigration model investigated in Ref. [16], where a new completely randomized species was introduced at each ''mutation step''. However, the present evolution model has substantially higher diversity in the high diversity state, provided that the parameters are chosen so that in the two models the interaction density c has the same value. The higher diversity in the evolution model is accompanied by many more small (order of size one) patches than the immigration model.
Finally, we would like to underline that the model implicitly incorporates an effective ''fitness'' with unusual properties. First of all, the fitness is clearly context dependent as the growth and the existence of a species entirely relies on the surrounding species. In this sense the model bears resemblance to the model of Ref. [27], where the stability of a species was defined in terms of its neighbors. Secondly, the exponential growth of a population in the current model is limited to a very short period, whereas a frozen steady state dominates most of the evolutionary trajectory of any species. We believe that our focus on short bursts and collapses captures the large scale evolution more accurately than a fitness defined by a potential for exponential growth.
The presented evolution scenario also suggests a reinterpretation of the Red Queen hypothesis of Van-Valen, who proposed that an observed constant species survival with geological time reflects a survival race where everybody has to improve just to maintain an unchanged survival chance [28]. We have observed that the survival of species lineages through time does not depend on obtaining large populations or on proliferation on a short timescale; instead, it depends on the extent to which a species hedges against hostile attack by splitting its population into isolated patches. Small populations of isolated species are often predicted to exist for very long time intervals. This is in contrast to the common observation that larger population indicates longer survival [29]. Of course, in reality, a species with too small population will not be stable; in the present model, one site is already assumed to be enough to sustain the species as long as there is no other species attacking it, indicating that one site is already above the minimum population size to avoid purely random extinction due to fluctuation of the population size. If the present mechanism of species survival is in effect, there can be a non-monotonic dependence of the species long-time survival on its population.
Although we have addressed the problem of non-motile species it is tempting to speculate about its overall behavior in the context of the fossil record of all animal species. It has been found in Ref. [30] that the long time survival of genera does not correlate with their short time ''success''. Although these data deal with duration of taxonomic orders with a short time ''success'' quantified by genera diversity, they also suggest a conceptual separation of the short time success from the long time survival.