Examining Community Stability in the Face of Mass Extinction in Communities of Digital Organisms

Paleontologists have long recognized that throughout the history of life on Earth, there has been a regular turnover of species (background extinction) punctuated by mass extinctions: infrequent but highly destructive crises that eliminate a sizable fraction of preexisting biological diversity, morphological variety, and ecological structure in geologically rapid episodes [10, 40, 74, 99]. These events are certainly remarkable for their destructive aspects; they are perhaps most notable for their marginalization or outright elimination of previously dominant incumbent taxa. Yet, mass extinctions are also a major creative evolutionary force. Removal of long-standing incumbents creates new opportunities for the subsequent evolution and diversification of surviving clades, often ones that were minor components of the pre-extinction biota. Processes of adaptive radiation and convergent evolution can then lead to new dominant taxa and new ecosystems [7, 72]. Mass extinctions can thus channel evolution in new directions unpredictable from previous background situations, permanently altering the ecological and taxonomic characteristics of biological communities [5, 53, 103].

In spite of the dramatic ecological and evolutionary consequences of mass extinctions and recoveries, they have been difficult to study because of problems posed by the uneven nature of the fossil record. Due to geological control of record quality, many previously studied paleobiologic patterns, including those of extinction and recovery, may be partly or wholly artefactual [43, 56]. Disappearances of taxa from the fossil record, due to geological hiatus or lack of suitable rock formations, may be incorrectly attributed to mass extinction until contrary evidence is discovered [7, 52, 101]. Different ways of adjusting data to compensate for incomplete preservation have resulted in differing conclusions about rates and processes involved in recovery from mass extinctions [4, 58, 66]. Beyond these problems, there is a lack of repeatability and control; one can never know how evolution might have proceeded under different conditions. For example, would a particular biological community have remained “stable” (by some measure) without a mass extinction and recovery? This last issue is worth considering in light of the phenomenon of coordinated stasis, where whole paleocommunities purportedly show taxonomic and ecological persistence over geologically significant time periods, interrupted by abrupt mass turnover [15–17, 21, 33, 34, 44, 48, 50, 77, 82, 87, 102]. More generally, we still do not know whether the consequences of different kinds of perturbations (instantaneous versus continuous, random versus selective) are different and repeatable. A deeper understanding of these issues is important not only for comprehending past mass extinction and recovery, but also because current global change represents a continuous, gradual anthropogenic perturbation [20, 73, 75].

Given the aforementioned problems with studying mass extinction and recovery from fossil data, it becomes highly desirable to have an experimental platform with which questions about these phenomena (and other outstanding general problems in long-term evolution) could be addressed without the complications of the geological record. Digital evolution offers some compelling advantages for studying community stability, and mass extinction and recovery more generally. First, and most importantly, it permits a truly manipulative, experimental approach to the problem. Its biggest strength is the ability to “replay life's tape,” to use the phrase of Stephen Jay Gould. One may use the capabilities of digital evolution to make changes at specific time points in what would otherwise be identical experimental replicates. For example, how would evolution unfold differently in the presence and absence of an extinction/recovery treatment? One may contrast paleobiological cases, which compare post-extinction and pre-extinction states (which are all one knows in the real world) with counterfactual cases, which compare the known post-extinction state against the corresponding state obtained with different treatments (e.g., no extinction). One may also examine how different kinds of model worlds respond to particular treatments, or vary the number and/or composition of survivors. Second, data may be recorded in far greater detail than is possible with the fossil record, and would be free of the preservational artefacts that complicate quantitative analysis of fossil data. Although not all the details and complexity of the geological record can be re-created with digital evolution, it is highly desirable to have a tractable system that allows one to focus on the key data sources that are most appropriate for investigating the problem at hand, and to design appropriate and tractable follow-ups. Third, because digital evolution is an actual instantiation of Darwinian evolutionary processes, it allows the problem to be addressed with greater realism of population-dynamic, ecological, and evolutionary factors than in other types of mass extinction simulations (e.g., [95]). As we are interested in both short-term and long-term effects of the perturbations we model (primarily the latter), our consideration goes beyond ecological scales and extends to paleontological time scales corresponding to thousands of post-extinction generations.

We use the digital evolution software Avida v.2.4.4 [2, 80] as our experimental platform, allowing a complete record of the course of evolution. An Avida world consists of a grid of digital organisms (hereafter referred to as Avidians), which in turn consist of a genome of instructions written in an assembly-like, Turing-complete computer language, along with a simple virtual computer chip that executes those instructions. The genome encodes the ability to self-replicate and carry out particular computations. Segments of code responsible for particular functions are analogous to genes. Avidians execute the programs encoded by their genomes, including commands that allow them to copy themselves and divide to produce a daughter organism. Memory for a daughter cell is allocated, and instructions are copied one at a time from parent to daughter. Copying is subject to point mutations, insertions, and deletions. These genomic mutations can indirectly affect an Avidian's ability to self-replicate or perform other computational functions, and their effects depend on their interaction with the rest of the genomic background. Most mutations are deleterious or neutral in their phenotypic effects, but a few are beneficial [61]. Thus, there is a genetic basis for adaptation and speciation. In our experiments, the Avidians are completely asexual, but this is no barrier to using them to address problems involving ecologically diversified populations, as asexual organisms can form phenotypically and functionally distinct phylogenetic clusters just as well as sexual organisms [6].

Each Avidian occupies a cell in the population grid, which is of finite size. Upon division, a daughter organism is placed in the population, killing any previous occupant of the destination cell (various placement options are available). This feature emulates removal of a healthy microbial cell from a chemostat, introducing an element of drift. All Avidians in a population receive a basal number of CPU cycles (analogous to energy), which enable their programs to run. However, Avidians can earn extra energy and accelerate their execution if they evolve the ability to perform certain computational functions. An Avidian that can perform any such functions receives an energetic reward by metabolizing a corresponding virtual resource into additional CPU cycles, enabling faster replication and creating variation among individuals in replication rate. Computationally, more complex functions are built by co-opting evolved mechanisms for simpler ones [62, 105]. Absolute time in Avida is measured in updates, during which an average of 30 instructions are executed per organism in the population. A generation typically requires 5–10 updates, with the precise number depending on the complexity of the Avidians present.

In our experiments, we used a 60 × 60 toroidal grid, initially seeded with a single copy of an asexual ancestor 50 instructions in length, capable only of self-reproduction. Mutations occurred during the copy procedure at a rate 0.0075 per instruction copied, and insertion/deletion mutations during division at a rate of 0.05 per divide, permitting change in both genotype and genome size.

The environment configuration used here contains two key features that enhance the model world's ecological realism. First, there are multiple depletable resources, in which use of a resource by one Avidian lowers the availability of that resource for other individuals requiring it. A low concentration of a resource reduces the benefit gained by performing the associated computation, favoring organisms that target underutilized resources. This resource competition model is analogous to that occurring in bacterial chemostats. Competition between Avidians that need exactly the same resources will favor organisms that use them most efficiently, but organisms with minimal overlap in resource requirements may coexist [23, 25]. Resources have no spatial structure; all organisms access the same resource pools regardless of location on the population lattice. Second, only a limited number of resources are supplied exogenously. The remainder occur as metabolic by-products when the Avidians consume available resources by successfully completing associated computations. These by-products then become resources to be metabolized by other Avidians, introducing a degree of ecological interdependence into the population dynamics (Figure 1). This second feature allows for evolution of rudimentary trophic structures through facilitative niche construction [79]. Although not a true food web because of the lack of true predation, evolution of such cross-feeding relationships has been observed previously in microbial evolution experiments [14, 57, 91]. The set of particular functions (and associated resources) used by any Avidian defines its ecotype, a concept analogous to an ecological guild: All Avidians that belong to the same ecotype perform the same functions and use the same resources, but do not necessarily descend from the same most recent common ancestor. We define ecotypes by the presence/absence of particular functions, but in practice apply the following screening procedure: If a particular function's execution makes up 5% or less of an Avidian's total functional output, it is set to zero, giving greater weight to those functions that are the major components of fitness for the Avidian.

In our experiments, we used an additive reward model. The energy gained for performing a rewarded function was added directly to the Avidian's current execution rate, and could be applied to multiple executions of that function over an organism's gestational cycle. Further information on how computational functions are evaluated and confer additional energy is provided in Supplement A in the online supplementary materials for this article at www.mitpressjournals.org/doi/suppl/10.1162/ARTL_a_00272.

Previous work [104, 106] investigated two principal types of mass extinction: pulse and press extinctions (Figure 2). Geologically, a pulse extinction happens with sufficient speed and power that adaptive change cannot occur during the extinction episode, although adaptation may take place afterwards. Conversely, a press extinction occurs over a longer period that allows for an adaptive response in affected populations [36]. In Avida, a pulse extinction is an instantaneous mass culling of individuals from the population, with survivors selected at random from the pool of all viable organisms. By contrast, a press extinction involves a press episode, a period of altered environmental conditions—in this case, greatly reduced inflows of basal resources—that persists long enough to allow an adaptive response. The aim is to trigger a massive drop in productivity, leading to bottom-up collapse of the ecological relationships in the community. Such a mechanism has been implicated in several major mass extinctions in the Earth's history [10, 40, 99]. Ecological recovery is then initiated by restoring resource inflows to pre-extinction levels. Unlike other models (e.g., [42, 97]) that used direct, targeted removal of individuals with particular phenotypic characteristics, here the ecosystem is allowed to adjust and evolve entirely on its own in response to the changed abiotic environment.

For the work described herein, the following experimental treatments were performed:

• 1. 

  Uninterrupted evolution (Control). Each replicate runs unperturbed for 205,000 updates.

• 2. 

  Strong Press episode. Each replicate runs for 100,000 updates before the press episode; these pre-extinction histories are identical to the Control. Resource inflows are then lowered by two orders of magnitude for 5,000 updates. This treatment is applied uniformly across experiments as absolute time, much as real extinction-driving crises act independently of biological generation times. Resource inflows are then restored for a subsequent 100,000 updates of evolution.

• 3. 

  Weak Press episode. Pre-extinction histories and press duration are as above, but the reduction in resources is not as severe. The low resource level was about tenfold higher than for Strong Press, so as to most closely match the phylogenetic attrition of Weak Pulse extinction ([104, 107]; see below).

• 4. 

  Strong Pulse extinction. Pre-extinction histories are as before. At 100,000 updates, an instantaneous mass cull of the population is performed. Survivors are picked randomly from among viable organisms, with no environmental alteration. This treatment is followed by 100,000 updates of recovery, again equal to the time for pre-extinction evolution. The survival rate of this strong cull is 4/3600 individuals (0.1%). This was the minimum number needed to ensure population continuity across all replicates; fewer survivors than this often resulted in population collapse due to stochastic effects or retention of highly unfit individuals.

• 5. 

  Weak Pulse extinction. Pre-extinction histories are as before, but the survival rate of the cull is 36/3600 individuals (1%). As Avidians may evolve to perform more than one computational function per organism, this is sufficient to preserve most of the community's ecological integrity.

We performed 100 replicates of each treatment. Each replicate with a particular random seed was identical in its evolutionary history up to 100,000 updates and diverged after that according to the treatment type.

A mass extinction and recovery may lead not only to a change in the total number of ecotypes present in the community, but possibly also to changes in community composition and structure. We address the following questions:

• 1. 

  To what extent were digital communities stable (by some measure) in the absence and presence of mass extinction?

• 2. 

  When ecotype identity was considered, did mass extinction and recovery cause greater rearrangement of dominant ecotypes than in corresponding Control replicates?

• 3. 

  Did mass extinction and recovery often result in the emergence of new stable states, a community with a persistent, differing taxonomic and ecological balance?

• 4. 

  When eco-phenotypic data suggested stability, was this supported once genealogical considerations were taken into account?

We use three different metrics with which to measure stability, each of which captures a different facet of the community. While they are not the ones most commonly used by ecologists, we chose them because of their suitability for the data generated by Avida and for ease of interpretability.

(a) Normalized mutual entropy (NME) as an index of specialization. Another concern in ecology, biogeography, and sociology is to what extent natural communities are generalized or specialized in some aspect, that is, the degree of division of labor of a set of individuals for a considered set of tasks (e.g., particular behaviors). For systems that may be clearly represented in the aforementioned manner, a promising approach is that of normalized mutual entropy (NME [38, 39]), the degree to which knowledge of an individual specifies knowledge of the tasks it performs. This is a somewhat more elaborate application of information-theoretic entropy-based measures to ecological data, a practice that has quite respectable precedent [68]. This metric has previously been applied in investigations of division of labor in social systems [12], but can be adapted for biogeographical studies [37]. It is particularly well suited to analyzing Avida data, since Avidian ecotypes in a population may be expressed as two-dimensional matrices, where each individual (rows) performs each rewarded function (columns) a certain number of times over the course of its life cycle. As there are more individuals than functions, the appropriate normalization is division of individuals into tasks (DIT), which is the mutual entropy of individuals and functions, divided by the marginal entropy of individuals [38]. In the formulation used here, the index is bounded between zero (no entropy, corresponding to a pure specialist community) and one (no information, corresponding to a pure generalist community). We use DIT as an overall index of specialization for the community, showing the extent of shifts in the degree of specialization over time. We can thus distinguish whether these shifts happen spontaneously due to internal community dynamics, or result from our experimental perturbations. As the calculation produces only a single point value of DIT for a given community, we generated nonparametric errors through repeated (5,000×) two-sample bootstrapping of the source community [32]. Details of the calculation of DIT, with a worked example, are shown in Supplement B in the online supplementary materials.

(b) Mean rank shift of the ten most common ecotypes (MRS₁₀). Ecologists have often sought to use abundance distributions for making mechanistic inferences about community structuring forces, but such approaches obscure species identity, and cannot detect shifts in dominance. Yet, such internal community dynamics are often the response of interest in manipulative experiments or long-term ecological monitoring, and metrics that incorporate information on compositional identity are highly desirable [70, 71]. An identity-sensitive metric used previously for assessing such community change is mean rank shift [24], which quantifies relative changes in species rank abundances over time. We use this metric to assess community change from different temporal perspectives, which we designate the prospective (from the past looking into the future) and retrospective (from the future looking into the past) comparisons. In the present context, the prospective comparison uses the immediate pre-treatment community as a reference point, while the retrospective comparison uses the end-experiment community state as the reference point. Since Avida communities contain many rare ecotypes, we modified this statistic to examine only the ten most common ecotypes in a community, which we designate as MRS₁₀. This limitation captures changes in the ecologically most important segment of the community, while filtering out most background noise. Further details of calculation, interpretation, and a worked example are given in online Supplement C. For comparison, we also use the more familiar Bray-Curtis (or Sorenson quantitative) index [68], which works with actual abundances rather than ranks; these results, which are completely congruent with those for MRS₁₀, are also presented in online Supplement C (Figure S5). We emphasize the MRS₁₀ results because they are more easily interpretable than the Bray-Curtis in terms of the number of ecotypes contributing to a shift of one standardized measurement unit (see the example in online Supplement C).

The DIT (including errors), MRS₁₀, and Bray-Curtis metrics were calculated with custom-written MATLAB scripts. Resampling of communities for DIT used functions from MATLAB's Parallel Computing Toolbox.

(c) Genealogical stability scores (GSS). Our DIT and MRS₁₀ data are derived purely from Avida's ecotypic output; there is an implicit assumption that for a particular ecotype at a given time, an ancestor-descendent connection exists between the most recent common ancestors (MRCAs) of that ecotype at both that time and earlier times. In order to verify whether or not there were enduring genealogical connections between ecotypes, we derived additional metrics that made use of genealogical information. The first metric is simply the proportion of shared ecotypes between the later and earlier times that had valid genealogical connections. The second is a more elaborate measure that incorporates information from ecotype, abundance, and genealogy. Briefly, when comparing two communities from different points in time that have substantial ecotype overlap and appear stable by one or both of MRS₁₀ or DIT, that stability is penalized if there is no genealogical connection between the MRCAs of an ecotype that is present at both times, that is, the ecotype at the later time is an ecological doppelganger of the one at the earlier time, but of different phylogenetic origin. The method, along with worked examples and an illustrative case, is described in detail in online Supplement D. For both these metrics, we compare the ecotypes at the end of the experiment with those in the immediate pre-treatment community (effectively a retrospective comparison).

Transformation of Avida population files into phylogenetic tree representations was performed with previously written software [107]. Additional labeling of tree files and filtering of Avida lineage and clade data were performed with custom-written MATLAB and R scripts. Tree visualization and MRCA identification were done with FigTree v.1.4.2 (http://tree.bio.ed.ac.uk/) and cross-verified with custom-written R scripts using functions from the ape [84], phytools [86], and Picante [55] packages.

Time series of ecological data may be used to infer periods of stability [51], either with reference to a preexisting level of the measured index [69], or by comparing alternative treatments for manipulative experiments [85]. We use both approaches here, particularly for DIT.

We identified stable periods within the MRS₁₀ and DIT time series using breakpoint analysis as implemented in the R package strucchange [108, 109]. This method has been employed previously in analysis of long-term ecological data to detect ecosystem regime shifts [67]. We use this procedure to investigate the effect of different external perturbations on change in community composition and division of labor. Comparison between different indices can reveal synchronicity and coordination of the changes at different levels, and evidence for stability. Where this method detected significant time series breaks, we always cross-checked the implied community changes against the ecotype lists for time points around the break (before and after), as the event underlying the break may be statistically but not ecologically significant.

We analyzed the DIT data as follows. First, we compared the pre-extinction and end-experiment average levels of DIT. The pre-extinction average level was determined by averaging the values from the immediate pre-treatment time to the first older breakpoint (see above) in the time series. Similarly, the end-experiment level averaged values from the end of the experiment to the first older breakpoint. To handle temporal autocorrelation, we first used iterative ARIMA modeling [26] to determine the optimal values for the autoregressive and moving-average parameters of the time series that would remove most of the autocorrelation. We then used a likelihood ratio test in order to compare the difference between the likelihoods of the ARIMA-adjusted data with those of a null model where the data were subtracted from the global mean of the pre-treatment and end-experiment time series. The aim was to see if the end-experiment state differed systematically from the pre-treatment state in a particular direction with respect to division of labor. We then repeated the above procedure, but this time for the end-experiment averages of the experimental treatments versus those of the corresponding Control replicate. Here, we wished to see if the end-experiment state of a treatment differed systematically from that of a Control. The iterative ARIMA procedure and likelihood ratio test were developed and implemented in R by T. Reitan, using R system library functions.

Second, we further developed a stability score based on examination of the entire post-treatment trajectory of each replicate, treating the Control trajectory as a reference level. We took this additional step because the previous analyses may obscure trends in DIT, particularly when a Control replicate itself is unstable. As this additional metric produced results completely congruent with the previously described analyses, we detail its calculation and results in Supplement B (Figure S2, Figure S3, Table S2) in the online supplementary materials.

We analyzed the MRS₁₀ data as follows:

• (a) 

  For prospective comparisons, we simply plotted the average time series for each treatment ± two standard errors; trends were sufficiently clear that formal statistical analyses were unnecessary. In order to test whether or not this metric remained static over time for Control, we employed Dunnett's test for multiple comparisons [31], using Control at 1,000 updates after treatment as a reference standard (we could not use the reference community itself, since all starting values for MRS₁₀ were zero); this test was performed for prospective comparisons only. Tests were performed in R using standard library functions; the Dunnett test used the glht( ) function from the Multcomp library [45].

• (b) 

  For retrospective comparisons, we noted considerable among-replicate variability in the results. Following treatment, some communities quickly settled into a state indicating low rank shift (and implied high similarity) compared to the end-experiment reference community, whereas others showed clearly more steplike traces when plotted. These results suggested that many populations recovering from extinction did not quickly form stable, mature post-extinction communities, requiring some way of assessing between-treatment variation in the number of intermediate community states. We therefore used recursive partitioning as implemented in the R system package rpart to fit step functions through the retrospective time series data. The plateaus in the predicted data were deemed locally stable stages (lss). If more than 95% of the observations for a particular replicate were less than the control end-experiment lower 95% confidence interval, that replicate was considered to have only a single state. To avoid overfitting of the data (and a resulting very large number of steps), all replicates were analyzed with control parameter settings of minpart = 40 and cp = 0.01.

All supporting data and R and MATLAB analysis scripts are available from Dryad (https://doi.org/10.5061/dryad.tq1g3j7). The source code for the versions of Avida and TreeLoader used to perform the experiments and convert population files into phylogenetic trees can be found at https://www.dropbox.com/home/Luo_etal_ALIFE_SourceCodes.

In this study, we used digital evolution as instantiated in Avida to investigate the stability of cross-feeding communities in response to mass extinction of two different types (Press and Pulse) and two different magnitudes (Strong and Weak). We investigated community stability using three different metrics based on division of labor (DIT), ecotype rank abundance (MRS₁₀), and genealogical connectedness (GSS). We found that:

• The analysis of DIT time series showed that the end-experiment state of Strong Press communities more often differed statistically and ecologically from both the corresponding pre-treatment and end-Control states than did the other treatments, but also had the widest variance of outcomes, with no systematic bias in either direction. Strong Press replicates also tended to establish a more stable post-treatment trajectory of DIT that differed decisively from the corresponding Control trajectory than did the other treatments. However, only Strong Pulse showed a statistically detectable bias towards communities that were more generalized than corresponding Controls. Control and Weak treatments showed few (if any) notable differences in these respects.

• Using MRS₁₀, we showed that a relatively small portion of the data set could be considered stable, with most (60–65%) Control replicates featuring between one and four replacements of dominant ecotypes over the length of the post-treatment period. Control showed a statistically significant signal of increasing MRS₁₀ over time, implying most replicates were not stable by this measure even without treatment. Strong Press communities had the largest and most consistent rank shifts of all treatments, Strong Pulse the second most, and Weak treatments generally did not differ significantly from Control.

• Both MRS₁₀ and DIT showed that recovery was more complicated from Strong treatments than from Control and Weak treatments. Strong treatments required more transitions in ecotype rank abundance for the final community state to take shape, whereas Control and Weak treatments were again quite similar in these characteristics, with fewer rank-abundance transitions and more total replicates that featured only a single locally stable stage (lss) over the length of the post-treatment period.

• Results using information from genealogy were generally consistent with those from analyses based only on eco-phenotypic data. For ecotypes that were shared between pre-treatment and end-experiment populations, Control and Weak treatments showed the greatest degree of genealogical continuity (regardless of rank or abundance), Strong Press the least, and Strong Pulse again intermediate between the two groups. When genealogical weightings were applied to a Bray-Curtis-like metric that included both ecotype and abundance information, Control and Weak treatments were all affected about equally, while the two Strong treatments showed greater reduction in stability, with Strong Press again most heavily affected.

Using DIT to investigate shifts in overall community specialization clearly illustrates differences in shifts due to internal dynamics versus those due to treatment-driven perturbations. Figure 3 shows two clear cases highlighting those differences. The community shown in panel (a) transitioned spontaneously from a generalized to a pure-specialist state before imposition of treatments, and in the Control replicate remained so for the duration of the simulation, indicating long-term stability of the specialist community. Of the treatments, only Strong Press disrupted this state, causing the emergence of a community dominated by more generalized organisms than either pre-treatment or end-Control states; DIT had still not settled into a new stationary state. All other treatments failed to disrupt the pre-treatment state, and the community remained dominated by specialists. In contrast, the community shown in panel (b) was more strongly generalized before treatments, and remained so in the absence of perturbation. This community was shifted into a more specialized state by the Strong treatments (particularly Strong Pulse), but recovered quickly to its previous state after the Weak treatments. Over the entire data set, using a criterion of a change of ≥0.1 in either direction to indicate ecologically meaningful change, we found among Controls 11 stable communities, 40 average communities, and 49 unstable communities.

The final state of DIT could be either more generalized or more specialized than the pre-treatment state (Figure 4). Overall, Control had the narrowest range and smallest variance of differences, while Strong Press had the largest (Table 1). However, only for Strong Pulse did the mean difference differ negatively and significantly from zero, indicating a slight tendency for Strong Pulse replicates to recover to a more generalized state than before treatment. The two Strong treatments had a wider variance of differences, differing significantly from Control, Weak treatments, and each other, while the Control and Weak treatments differed only from the Strong treatments. These same tendencies were seen when comparing the end-treatment states with the end-Control state (Supplement B, Table S1; Figure S1 in the online supplementary materials), and when considering the entire post-treatment time series (online Supplement B, Figure S2, Figure S3, Table S2). Except for Strong Pulse, the other treatments showed no systematic bias in end-treatment DIT state versus either pre-treatment or end-Control states.

Using a change of ≥0.1 (in either direction) as the criterion for ecologically meaningful change, we identified among the Controls 11 stable communities, 40 average communities, and 49 unstable communities. There were no significant differences between pre- and post-extinction states for any treatments in stable communities. For average communities, only some replicates yielded significant differences, and differences were always significant in unstable communities.

As mentioned previously, we used a criterion of 0.1 as a threshold for stability. We could thus divide the results into communities that were stable (almost always around 0.1 or less), average (between 0.1 and 0.3, up to four replacements of dominant ecotypes), or unstable (0.3 or greater, four or more replacements of dominant ecotypes). Summary results are in Table 2 and Figure 5.

There was significant variation among treatments, with Strong Press in particular having an excess of average and unstable communities (from both perspectives) compared to Control; at the other extreme, Weak Pulse hardly differed at all. The effect of Weak Press was equivocal; while there were fewer stable communities and somewhat more average and unstable communities than for Control, it was not nearly as dramatic as for the Strong treatments. The response of Strong Pulse was intermediate between Strong Press and the Weak treatments, as is more easily seen in the prospective view.

Clearly significant variation among treatments in MRS₁₀ was shown for prospective comparisons over nearly all of the post-treatment period (Figure 5). Control communities close in time to the beginning of the post-treatment period did not differ from the (near-)reference community, but on average became increasingly dissimilar over time (Table 3).

Comparing between communities, the largest differences were produced by Strong Press, which differed strongly from all other treatments (Figure 5, dark blue traces). The Control and Weak treatments all showed comparable values and did not differ significantly from each other over the post-treatment interval (except for a brief window very early in the recovery for Weak Press). On average, Weak Press was around the upper 95% confidence intervals of Control and Weak Pulse (Figure 5a, cyan traces), and it had fewer stable replicates than did the former two treatments (Table 2). Strong Pulse showed an intermediate response, but remained statistically distinct from all other treatments (Figure 5a, red traces). These results show that, on average, the Strong treatments resulted in an intermediate to large perturbation of pre-extinction community structure. The Control and Weak treatments resulted in minimal perturbation of community structure (somewhat stronger in the case of Weak Press), with dominant ecotypes generally remaining so. However, even in Control replicates, the large number of average communities (more than half of all replicates: Table 2) and steadily increasing average MRS₁₀ over time (Table 3) imply an overall lack of community stability even for this (non-)treatment.

As mentioned previously, for retrospective comparisons, we found a considerable diversity of outcomes for changes in individual replicates, even in Controls (Figure 6). Some cases were quite stable, maintaining low values over the entire post-treatment time period (Figure 6a). In such communities, dominant ecotypes generally maintained their pre-treatment relative ranks, with very little introgression of new types. However, there were also unstable cases, where substantial rank shifts occurred quickly and unpredictably (Figure 6b). Where large rank shifts occurred, at least two pre-treatment dominant ecotypes were replaced by ones with functional combinations that were previously rare or absent; many such communities were mosaics of old and new ecotypes, but complete turnover did not occur. Responses for Strong treatments, however, ranged from intermediate turnover to complete renewal of the dominant types. For prospective comparisons, many communities quickly stabilized at higher MRS₁₀ values, implying that the post-extinction community state differed substantially and permanently from the pre-treatment state (Figure 5, dark blue trace; Figure 6c, d, stippled traces). For the corresponding retrospective comparisons, rank shifts could soon stabilize or occur in a more stepwise manner over the post-treatment period (e.g., Figure 6c, solid trace). Our interpretation of such results is that the formerly dominant pre-treatment ecotypes either quickly went extinct or became very rare. By contrast, the composition of the post-treatment community could take shape either rapidly or in stages over the course of the recovery period, with the largest shifts driven by multiple new dominant types that were ecologically marginal or absent in the pre-treatment community. We analyzed this more closely through determining lss for retrospective comparisons (Table 4).

We found significant variation among treatments for the number of lss of retrospective MRS₁₀ (F_(4,495) = 19.04, p < 1.41 × 10⁻¹⁴). The results were consistent with those for prospective comparisons; there were no significant differences in lss among Control and Weak treatments, Strong Press had on average the most lss, and Strong Pulse was again intermediate between the two groups. When considered on a within-treatment basis, a larger fraction of the data for Control and the Weak treatments showed only a single lss, Strong Press had the fewest such replicates, and Strong Pulse was again intermediate. Seen this way, the Control and Weak treatments tended during recovery to most quickly establish a community state that resembled the end-experiment one (if the immediate pre-treatment community was at all perturbed). Strong Press communities went through the most transitional states before establishing a final state (indeed, many were still in transition at the end of the experiment), while Strong Pulse fell between these extremes.

The mass extinction treatments affected not only ecological stability as indicated by the preceding metrics, but also the underlying phylogenetic and genealogical structure of the affected populations. The first, simple measure of genealogical stability revealed that markedly fewer ecotypes shared between the end-experiment community and the immediate pre-treatment community maintained genealogical continuity across the treatment and recovery period for the two Strong treatments than for the other treatments (Figure 7a). The fraction of shared ecotypes with connections was quite similar for Control and Weak treatments, with about 78%, 75%, and 72% for Control, Weak Pulse, and Weak Press respectively. This fraction fell to about 53% for Strong Pulse and 42% for Strong Press, indicating an increased occurrence in these treatments of shared ecotypes that were convergently evolved ecological doppelgangers of their pre-treatment counterparts. Importantly, the differences among treatments were in the same direction as MRS₁₀ and statistically significant in the same way (one way ANOVA, F_(4,495) = 98.36, p < 1.5 × 10⁻⁶¹), with significant differences between the Strong treatments after multiple comparison testing.

The second metric, which incorporated information from ecotype identity, abundance, and genealogy, demonstrated that including genealogical information could alter the value of the metric—and thus inference of the degree of community stability—versus the assumption of connections between all ecotypes (Figure 7b). When unweighted, essentially the same conclusions as using MRS₁₀ (Figure 5) or the Bray-Curtis index (online Supplement C, Figure S5) obtain, but these are changes due only to ecotype turnover. Including genealogical weightings reduced the value of the metric further still for all treatments, by an average of 7.3% to 8.2% for Control and Weak treatments, about 11% for Strong Pulse, and about 14% for Strong Press (Figure 7b).

Our results provide information about a general pattern in how these cross-feeding digital communities are structured. The generally small values of MRS₁₀ (Table 2) and stable values of DIT (Table 1) observed in the Control treatment suggest that most replicate populations contain an ecological core: a set of organisms that provides most ecosystem services and ecological structure for the population. Around these core organisms exist many less common, often transient types generated by mutation. This is analogous to natural ecosystems that feature resident and tourist species, the former characterized by long-standing persistence and high abundance, the latter by infrequent occurrence and low abundance [70], a key difference being that in real communities on ecological time scales, nearly all tourists enter the community through migration, rather than in situ through speciation. Large shifts in either MRS₁₀ or DIT indicate an alteration of that core, usually through replacement of several residents by former tourist types. Such shifts were observed mainly in the Strong treatment communities. However, even undisturbed Control communities could also often experience substantial ecotype turnover (e.g., Figure 6b), as demonstrated by the significant signal of increasing departure (Table 3) from near-pre-treatment rank abundance (though not to the extent of the Strong treatments). This may result from delayed initial evolution of ecotypes with top-tier functions (particularly the computationally more difficult xor and equals functions), or replacement of once-dominant ecotypes by new ones with altered functional groupings, even after all functions have evolved. One or more tourist ecotypes may displace previous core members, though it may take several independent appearances of these tourists for that to occur, as the initially rare types must first dodge drift [60], and be at least equal (or superior) in fitness to one of the existing dominant types. It is unclear if this slower turnover of ecotypes is due more to Red Queen phenomena [63] or community drift processes [46, 64]. We hypothesize that such replacements might alter relative fitnesses and ecological relationships of other community members, rearranging prevailing ecological structure and influencing subsequent evolution of the community [3, 8], and we plan to test this directly in future work.

On average, the Weak extinction treatments were insufficiently disruptive to produce large changes beyond those of Control populations, though their effect was somewhat more substantial in the case of Weak Press (Figure 5). In stark contrast, the Strong treatments upset the previous ecological core, allowing for the evolution and establishment of a new set of core types. This was most clearly evident in the case of Strong Press: It had the highest average prospective MRS₁₀ (Figure 5, Table 2), its communities required on average the most transitions between locally stable post-treatment states and generally passed through the most intermediate steps before settling into a terminal-like state (Table 4), it produced the largest range of differences in DIT (Figure 4, Table 1), and it had the most reduced genealogical stability (Figure 7). This treatment, as we have implemented it, acts as more than a simple filter of organisms, as do the Pulse treatments: It often resulted in near-total ecosystem collapse, effectively erasing the previous core at functional-genomic, ecological, and phylogenetic levels [105–107]. Ecosystem reestablishment started nearly from scratch, from a phenotypically altered, functionally reduced pool of Avidians better adapted to the harsh conditions of the Strong Press episode. This difference in post-treatment evolutionary landscape is also reflected in the relative contributions of genealogically unconnected ecotypes in Press versus Pulse (Supplement D, Table S4 in the online supplementary materials). The greater contribution from unconnected types in Pulse treatments demonstrates that, even with near-total annihilation of the pre-treatment community, the unaltered survivors in many cases retain the genetic and functional potential to regenerate lost ecotypes with the same ecological profiles, though there is also of course room for vacated ecospace to be occupied by newly evolved organisms with different ecological profiles. The kinds of changes resulting from Press treatments, on the other hand, favor replacement of marginalized or lost ecotypes by functionally different ones (which may be ecologically similar, but not identical), leading to a more thoroughly reshaped and restructured community.

In contrast to the highly visible differences created by Strong Press, Strong Pulse was generally intermediate in most aspects between Strong Press and Control on one hand and the Weak treatments on the other (Figure 5, Tables 1, 2, 4). In our implementation of pulse extinction, survival is random, and a population may be reestablished by a subset of Avidians ranging anywhere from all common, core organisms at one extreme to all rare, tourist organisms at the other. Previous work [104] showed a probability of 0.49 that at least two of the Strong Pulse survivors would come from among the two most common ecotypes. Further, the overall generalized nature of most replicate populations (average pre-treatment DIT = 0.465 ± 0.011 (s.e.), about 2–3 functions per Avidian) means that even a sample as small as those used for Strong Pulse may actually still contain many of the functions needed to rebuild the ecology (albeit in heavily reduced amounts), free of the eroding effects of a press episode. Survivors belonging to previously dominant groups may then quickly re-assume key ecological roles. Thus, communities recovering even from a Strong Pulse extinction might retain some “memory” of their previous state.

This is supported by the lower prospective MRS₁₀ than for Strong Press (Figure 5), and, more interestingly, the fact that among the distributions of DIT difference scores (Figure 4; Table 1), only the mean of Strong Pulse was significantly different (albeit small) from zero, indicating a slight tendency towards more generalized post-treatment communities—a result shown by no other treatment. For Control and Weak treatments, the numerical differences are simply not large enough (even if statistically significant: Figures 4, 5; Tables 1, 2, 4) to imply substantial regular ecological differences between the communities. Under Strong Press, destruction of the pre-treatment community is often so thorough that there is effectively no relation between the pre- and post-treatment communities with respect to DIT. The post-treatment community is just as likely to be more specialized as more generalized, as shown by the wider variance in DIT comparisons (Table 1; Figure 4), as well as the signed stable trajectory score (Supplement B, Table S2 in the online supplementary materials). These results highlight an important difference between the random Pulse treatment and the selective Press treatment. Generalist survivors may well have an early advantage in post-treatment expansion, as they can use multiple resources and may in many cases be able to survive by self-feeding. However, the difference here seems to be more than just a simple short-term expansion of lucky generalized survivors, in that Strong Pulse actually leaves a small, but enduring, signature in the community's functional makeup.

We were also concerned with whether our digital communities shifted into new stable states following a large disruption, that is to say, a community characterized by different ecological and taxonomic makeup, even with the same underlying abiotic environment, that does not easily revert to the previous one [29, 49, 85, 93]. For this digital system, the possibility of post-extinction alternative states was raised in previous work [106], based on simple phenotypic diversity-through-time results. We have explored it more fully here, as raw phenotypic diversity can indicate recovery to similar diversity levels but hide significant changes in community makeup. Our results show that sufficiently strong disruption could indeed often drive these digital communities into new states with a different taxonomic and ecological balance, shown by high prospective MRS₁₀ indicating a heavily altered set of core ecotypes, substantially altered division of labor as revealed through DIT, or both.

However, there are two key caveats to these results. The first is that while strong mass extinction often succeeded in driving communities into new states, we cannot say that they consistently represented alternative stable states—the assumption being that the immediate pre-treatment community itself represented a stable state. Indeed, our data indicate that stable pre-treatment communities were relatively uncommon. The overall signal of increasing prospective MRS₁₀ for Control (Table 3), as well as the substantial number of Control communities that were deemed average and unstable (Table 2) or showed a significant difference in DIT from the pre-treatment state (Table 1), demonstrate that even without perturbation many communities were still subject to change in the second half of the experiment. The point where the treatment was imposed divided the experiment into convenient before and after phases, but the community itself need not have been stable at that point (e.g., Figure 6b). Hence, it is more accurate to say that the treatments (especially the Strong ones) drove many communities into a state that differed substantially from that achieved by the Control, regardless of Control stability.

The second consideration is that even when new stable states emerged, they often did not result quickly following a Strong treatment. The higher average number of lss in retrospective MRS₁₀ and the low number of replicates that quickly converged to a community state with long-term stability (Table 4) indicate that after a Strong treatment, many early post-treatment communities are unstable and not yet mature. A persistent, stable community state may only emerge after a prolonged period of time, after passing through several transitional states (e.g., Figure 6c). Such early post-treatment instability was implied by changes in phylogenetic properties of these communities [107], and we have now demonstrated that more directly. Not all post-treatment instability resulted in transition to an alternative state: There were certainly a number of communities where early instability was followed by an eventual return to pre-treatment Control values. However, caution must be exercised with this interpretation, particularly depending on what metric is used to infer stability. Different types of communities may produce similar DIT values but have different ecotype composition. For any given replicate experiment, it is better to infer stability from consistency across different metrics—for example, if a return to pre-treatment values of DIT is accompanied by low rank shift, also indicating return to pre-treatment ecotype composition and ranking. Explicit consideration of phylogeny/genealogy then becomes necessary to determine whether ecotypes shared between pre- and post-treatment communities really are genealogically connected.

Since we designated ecotypes based on the presence or absence of rewarded computational functions, the possibility existed that new Avidians with the same ecotype as those from an already existing clade, but of different phylogenetic origin, could arise and replace their eco-equivalents during the course of an experiment, an occurrence that would go undetected. Additionally, the computational functions of Avidians can be re-evolved in lineages where they once existed and were lost due to the erosive effects of the press episode [105]. Hence, we had to establish whether there were indeed genealogical connections between the same ecotypes at widely separated points in time.

Both our measures of genealogical stability revealed that in all treatments, changes in genealogical connections further impacted stability beyond ecotype turnover alone (Figure 7). For all treatments, the genealogical changes alone do not contradict the ecological results (which would make the latter meaningless); neither do they cancel the ecological results, since they occur in the same direction and compound rather than counteract ecological changes. Although Control and Weak treatments differed little from each other, inclusion of genealogical information revealed that while Weak Press had about the same fraction of shared ecotypes as did Weak Pulse, its difference from the unweighted metric was slightly smaller, indicating less of an impact from genealogical changes than for Weak Pulse. Close inspection of the data revealed that in Weak Press, a smaller fraction of total replicates had at least one genealogically unconnected shared ecotype as compared to in Weak Pulse. Instead, there was a somewhat greater fraction of replicates where new dominant ecotypes arose from extreme rarity (Supplement D, Table S4 in the online supplementary materials). In Weak Pulse, on the other hand, a slightly larger fraction of replicates had at least one genealogically unconnected shared ecotype, and fewer with new ecotypes arising from extreme rarity (online Supplement D, Table S4).

Taken together, these results suggest that the Press and Pulse treatments exert their effects in different ways. The (slightly) more ecologically extensive changes seen in Weak Press are due more to press-driven ecotype turnover (and associated ecological restructuring) that has a lesser phylogenetic impact on shared ecotypes (but see [107]). By contrast, in Weak Pulse the phylogenetic impact on shared ecotypes is somewhat greater, with more convergent re-evolution of lost ecotypes and resulting lower ecological impact (which usually hardly differs from that in Control). The aforementioned tendencies also occurred in the Strong treatments, as both MRS₁₀ and unweighted GSS showed that the major effect of Strong Press is due to extensive ecotype turnover. Including genealogical information shows that the radical changes seen in most Strong Press communities are in fact due to the combined effect of both this ecotype turnover and tremendous phylogenetic attrition [107]. While the ecological restructuring wrought by Strong Pulse is generally less than that by Strong Press, the additional contribution from genealogy is nearly as large, also consistent with previous results that examined only phylogeny [107]. Hence, for a system such as this one where ecotype identity is defined through performance of particular ecological functions, additional consideration of genealogical/phylogenetic information can further illuminate the degree to which apparent ecological stability conceals an underlying phylogenetic dynamism.

One of our original goals was to determine if these digital communities showed sufficient evidence for stability that would result in coordinated stasis. While our digital data differ in a number of ways from available paleontological data, we can still assess coordinated stasis as a pattern through relative stability of rank abundance, community specialization, and genealogical connectedness among ecotypes over appreciable spans of simulation time. When only our eco-phenotypic metrics are considered together, only a fairly small fraction of Control replicates (25%) could be considered to show coordinated stasis, inferred through a combination of very low MRS₁₀ and near-constancy of DIT across the whole second half of the experiment, but this becomes less convincing after allowing for genealogical connectedness. Thus, we must conclude that the experimental conditions used here lack the necessary ecological complexity to consistently produce sufficiently stable communities where even the pattern of coordinated stasis is prevalent. Rather, our data better support an alternative paradigm of structural continuity [76, 88, 89] for these communities: Community structure is more or less persistent, but the particular actors (species, ecotypes, etc.) that fill particular ecological roles come and go dynamically. Our results, and the conclusion of the prevalence of structural continuity, are also consistent with previous results using the same digital system, with a different experimental and analytical setup [81]. However, given the existence of capabilities for movement and spatial heterogeneity of resources in Avida [35], examining the current leading hypothesis for coordinated stasis—habitat tracking [15, 18]—is an interesting and plausible direction for future work.

Previous work has covered the broad qualitative similarities of the results produced by this digital system to major extinction and recovery events from Earth history [106]. While many of the particulars of our current results are likely to be system-specific, we can draw further broad parallels based on the present results and advances that have been made in paleontology since that time.

• Initially immature post-extinction communities and stepwise assembly. Our results showed that rapid formation of a new community that was stable over the whole post-treatment interval was relatively uncommon for the Strong treatments; a more stepwise dynamic, indicating several different intermediate stages before emergence of the end-experiment community, was more common. Real post-extinction paleocommunities can also display such a variety of dynamics, as the component species and their ecological relationships change and stabilize in a succession-like process [47]. Even if overall trophic structure recovers quickly [94, 100], communities from the earliest stages of recovery may be atypical in comparison with later, more persistent community states, again best observed in post-Permian Triassic communities in both terrestrial [9, 90, 92] and marine [94] contexts, but also evident in the early Paleogene from a number of independent lines of evidence (reviewed in [30]). However, clearly stepwise recovery may depend on geographic context and may not be general at a global scale (compare [11, 22] versus [94]). In our digital system, these post-extinction turnovers and transitions may be due to delayed re-evolution of new functions [106] or a superior new implementation of already existing functions. The Cretaceous terrestrial revolution and Early Eocene Climate Optimum (reviewed in [47]) provide real paleontological analogues.

• Radical community turnover after devastating extinction, and new stable states. Our results also showed that Strong Press extinctions in particular resulted in the greatest tendency to produce new ecosystems that differed most substantially in ecological and genealogical makeup from both the pre-treatment and end-Control states, with substantial contribution from ecotypes that were previously marginal or absent. Such extensive restructuring of communities is a hallmark of the most ecologically disruptive mass extinctions and recoveries of the Phanerozoic, particularly the end-Permian [28, 101] and end-Cretaceous [1, 30] extinctions (although the latter is usually interpreted as a canonical pulse extinction [27]), which led to phyletically different sets of taxa becoming dominant members of the biota and the emergence of new ecological structures. The less-extensive restructuring produced by Strong Pulse extinctions more often resulted in post-treatment communities made up of a variety of pre-treatment survivors, genealogically unconnected re-evolved ecotypes, and formerly rare or new ecotypes with increased ecological prominence. The arrival and establishment of these new ecological players would act to prevent a return to previous states [83, 98]. The early Triassic [96] and early Paleogene [30] recoveries provide perhaps the best analogies, since while many of the available niches were filled through diversification from surviving taxa, others continued to be occupied by carryover taxa (e.g., Lystrosaurus across the Permo-Triassic boundary), or were filled by radiation from lineages that were related to, but had previously been less prominent than, the previous occupants—for example, the replacement of ecologically diversified archaic bird lineages by the modern neornithine clade in the early Paleogene [65]. The Weak treatments (particularly Weak Pulse), which generally did not produce significant departure from Control, are perhaps more reminiscent of the end-Ordovician, which, while devastating to organisms on lower taxonomic levels, did not eliminate many higher taxa (order or more inclusive) and did not have the well-documented long-term impacts of the end-Permian and end-Cretaceous events [41].

• Shift to increased generalization in Strong Pulse communities. Our DIT results indicated a tendency for a slight (but still detectable) shift to more generalized communities after Strong Pulse. This result has some interesting analogies with the paleontological record, where certain mass extinctions (end-Permian [92], end-Triassic [19], and end-Cretaceous [59]) have resulted in post-extinction communities with increased ecological generalization and/or biogeographic cosmopolitanism. Such a shift is projected for the modern biota, given ongoing anthropogenic extinction [13, 54, 78]. Unlike these paleontological examples, in which the increased generalization/cosmopolitanism eventually declines, our data show a (marginally) persistent effect of the extinction and recovery on division of labor.

Our results illustrate the potential for digital evolution to make many aspects of studying mass extinction, recovery, and community stability more transparent than is possible with the fossil record, and treat them in a truly experimental manner. We can study how evolution would proceed under different sets of circumstances and connect observed differences to the types of manipulations. We have shown that different metrics measure different facets of community change, and demonstrated that applying several approaches yields a fuller understanding of community change across multiple levels. This gives us confidence that further use and development of evolving digital systems, combined with measurement approaches borrowed from community ecology, can yield additional insight into problems in paleobiology and long-term evolution that are not amenable to other experimental approaches.

While one must exercise great caution in extrapolating results from a digital system like ours to the real world, we believe that our results have bearing on potential future trends for the ecology and evolution of the world's biota. It has become increasingly accepted that we are in the midst of an anthropogenic sixth mass extinction [20, 73, 75]. Human activities driving the extinction can be considered a continued environmental pressure that favors the survival of certain taxa while threatening and ultimately removing others, which our Press extinctions are analogous to. A key implication of our study is that full recovery from even a relatively mild biotic crisis often entails a transformed community at both ecological and phylogenetic levels, the latter being more pronounced with increasing ecological severity. If the biotic world comes to be dominated by taxa that both survive and thrive in conditions created by human activity, then, depending on the severity of the effects (whether weak or strong), recovery of the biosphere may ultimately result in radically altered ecosystems highly divergent from those we know today, with new taxonomic and ecological balances—as has already happened after the most severe mass extinctions in Earth history—that are difficult (if not impossible) to reverse and that, depending on their component species, may not be helpful to human well-being.

We thank two anonymous reviewers for comments and feedback that helped us improve the manuscript. This study was supported by research grants #31470435 from the National Natural Science Foundation of China (http://www.nsfc.gov.cn) and #KYRC201301 from the Chinese Ministry of Education (http://www.moe.edu.cn) to G.Y. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. G.Y. gives special thanks to Carlos Alberto Arnillas (Department of Biology, University of Toronto–Scarborough, Toronto, Canada) for suggesting the stable trajectory approach to the mutual entropy data.