Adaptive Phenotypic Plasticity Stabilizes Evolution in Fluctuating Environments

Fluctuating environmental conditions are ubiquitous in natural systems, and populations have evolved various strategies to cope with such fluctuations. The particular mechanisms that evolve profoundly influence subsequent evolutionary dynamics. One such mechanism is phenotypic plasticity, which is the ability of a single genotype to produce alternate phenotypes in an environmentally dependent context. Here, we use digital organisms (self-replicating computer programs) to investigate how adaptive phenotypic plasticity alters evolutionary dynamics and influences evolutionary outcomes in cyclically changing environments. Specifically, we examined the evolutionary histories of both plastic populations and non-plastic populations to ask: (1) Does adaptive plasticity promote or constrain evolutionary change? (2) Are plastic populations better able to evolve and then maintain novel traits? And (3), how does adaptive plasticity affect the potential for maladaptive alleles to accumulate in evolving genomes? We find that populations with adaptive phenotypic plasticity undergo less evolutionary change than non-plastic populations, which must rely on genetic variation from de novo mutations to continuously readapt to environmental fluctuations. Indeed, the non-plastic populations undergo more frequent selective sweeps and accumulate many more genetic changes. We find that the repeated selective sweeps in non-plastic populations drive the loss of beneficial traits and accumulation of maladaptive alleles, whereas phenotypic plasticity can stabilize populations against environmental fluctuations. This stabilization allows plastic populations to more easily retain novel adaptive traits than their non-plastic counterparts. In general, the evolution of adaptive phenotypic plasticity shifted evolutionary dynamics to be more similar to that of populations evolving in a static environment than to non-plastic populations evolving in an identical fluctuating environment. All natural environments subject populations to some form of change; our findings suggest that the stabilizing effect of phenotypic plasticity plays an important role in subsequent adaptive evolution.

1. Introduction

Natural organisms employ a wide range of evolved strategies for coping with environmental change, such as periodic migration (Winger et al., 2019), bet-hedging (Beaumont et al., 2009), adaptive tracking (Barrett and Schluter, 2008), and phenotypic plasticity (Ghalambor et al., 2007). The particular mechanisms that evolve in response to fluctuating environments will also shift the course of subsequent evolution (Wennersten and Forsman, 2012; Schaum and Collins, 2014). As such, if we are to understand or predict evolutionary outcomes, we must be able to identify which mechanisms are most likely to evolve and what constraints and opportunities they impart on subsequent evolution.

In this work, we focus on phenotypic plasticity, which can be defined as the capacity for a single genotype to alter phenotypic expression in response to a change in its environment (West-Eberhard, 2003). Phenotypic plasticity is controlled by genes whose expression is coupled to one or more environmental signals, which may be either biotic or abiotic. For example, the sex ratio of the crustacean Gammarus duebeni is modulated by changes in photoperiod and temperature (Dunn et al., 2005), and the reproductive output of some invertebrate species is heightened when infected with parasites to compensate for offspring loss (Chadwick and Little, 2005).

Evolutionary biologists have long been interested in how evolutionary change is influenced by phenotypic plasticity because of its role in generating phenotypic variance (Gibert et al., 2019). The effects of phenotypic plasticity on adaptive evolution have been disputed, as few studies have been able to observe both the initial patterns of plasticity and the subsequent divergence of traits in natural populations (Ghalambor et al., 2007, 2015; Wund, 2012; Forsman, 2015; Hendry, 2016). In changing environments, adaptive phenotypic plasticity provides a mechanism for organisms to regulate trait expression within their lifetime, which can stabilize populations through those changes (Gibert et al., 2019). In this context, the stabilizing effect of adaptive plasticity has been hypothesized to constrain the rate of adaptive evolution (Gupta and Lewontin, 1982; Ancel, 2000; Huey et al., 2003; Price et al., 2003; Paenke et al., 2007). That is, directional selection may be weak if environmentally-induced phenotypes are close to the optimum; as such, adaptively plastic populations may evolve slowly (relative to non-plastic populations) unless there is a substantial fitness cost to plasticity.

Phenotypic plasticity allows for the accumulation of genetic variation in genomic regions that are unexpressed under current environmental conditions. Such cryptic (“hidden”) genetic variation can serve as a source of diversity in the population, upon which selection can act when the environment changes (Schlichting, 2008; Levis and Pfennig, 2016). It remains unclear to what extent and under what circumstances this cryptic variation caches adaptive potential or merely accumulates deleterious alleles (Gibson and Dworkin, 2004; Paaby and Rockman, 2014; Zheng et al., 2019).

The “genes as followers” hypothesis (also known as the “plasticity first” hypothesis) predicts that phenotypic plasticity may facilitate adaptive evolutionary change by producing variants with enhanced fitness under stressful or novel conditions (West-Eberhard, 2003; Schwander and Leimar, 2011; Levis and Pfennig, 2016). Environmentally-induced trait changes can be refined through selection over time (i.e., genetic accommodation). Further, selection may drive plastic phenotypes to lose their environmental dependence over time in a process known as genetic assimilation (West-Eberhard, 2005; Pigliucci, 2006; Crispo, 2007; Schlichting and Wund, 2014; Levis and Pfennig, 2016). In this way, environmentally-induced phenotypic changes can precede an evolutionary response.

Phenotypic plasticity may also “rescue” populations from extinction under changing environmental conditions by buffering populations against novel stressors. This buffer promotes stability and persistence and grants populations time to further adapt to rapidly changing environmental conditions (West-Eberhard, 2003; Chevin and Lande, 2010).

Disparate predictions about how phenotypic plasticity may shift the course of subsequent evolution are not necessarily mutually exclusive. Genetic and environmental contexts determine if, and to what extent, phenotypic plasticity promotes or constrains subsequent evolution. Figure 1 overviews how we might expect different forms of phenotypic plasticity to result in different evolutionary responses after an environmental change. In Figure 1A, we would expect non-plastic populations to experience strong directional selection toward the new optimum (O_E2) after the environment changes. We would expect an adaptively plastic population (Figure 1B) to remain relatively stable after the environment changes, as plasticity shifts organisms' phenotypes to the new optimum. In Figure 1C, we would expect the non-adaptively plastic population to experience strong directional selection on their response to the new environmental conditions; indeed, such maladaptive plasticity may even put the population at risk of extinction in the absence of beneficial mutations.

FIGURE 1

Figure 1. Hypothetical reaction norms for populations comprising genotypes placed in different environments. A reaction norm describes phenotypic change (or lack thereof) induced by environmental variation (West-Eberhard, 2008). In all panels, two environmental conditions (denoted E₁ and E₂) are shown on the x-axis. The y-axis indicates the phenotype expressed in each environment with O_E1 and O_E2 designating the optimal phenotype for E₁ and E₂, respectively. Each pair of points connected by a solid black line denotes a genotype, with the points themselves representing its hypothetical phenotypes in each environment. We present three scenarios for how populations could respond to a change from E₁ to E₂. (A) A non-plastic population where phenotypes do not change with environmental shifts. (B) An adaptively plastic population where phenotypes dynamically adjust to the new optimum. (C) A population exhibiting non-adaptive plasticity where environmental change induces phenotypes further away from the optimum.

Experimental studies investigating the relationship between phenotypic plasticity and evolutionary outcomes can be challenging to conduct in natural systems. Such experiments would require the ability to irreversibly toggle plasticity followed by long periods of evolution during which detailed phenotypic data would need to be collected. Digital evolution experiments have emerged as a powerful research framework from which evolution can be studied. In digital evolution, self-replicating computer programs (digital organisms) compete for resources, mutate, and evolve following Darwinian dynamics (Wilke and Adami, 2002). Digital evolution studies balance the speed and transparency of mathematical and computational simulations with the open-ended realism of laboratory experiments. Modern computers allow us to observe many generations of digital evolution at tractable time scales; thousands of generations can take mere minutes as opposed to months, years, or millennia. Digital evolution systems also allow for perfect, non-invasive data tracking. Such transparency permits the tracking of complete evolutionary histories within an experiment, which circumvents the historical problem of drawing evolutionary inferences using incomplete records (from frozen samples or fossils) and extant genetic sequences. Additionally, digital evolution systems allow for experimental manipulations and analyses that go beyond what is possible in wet-lab experiments. Such analyses have included exhaustive knockouts of every site in a genome to identify the functionality of each (Lenski et al., 2003), comprehensive characterization of local mutational landscapes (Lenski et al., 1999; Canino-Koning et al., 2019), and the real-time reversion of all deleterious mutations as they occur to isolate their long-term effects on evolutionary outcomes (Covert et al., 2013). Furthermore, digital evolution studies allow us to directly toggle the possibility for adaptive plastic responses to evolve, which enables us to empirically test hypotheses that were previously relegated to theoretical analyses.

In this study, we conduct digital evolution experiments to investigate how the evolution of adaptive phenotypic plasticity shifts the course of evolution in a cyclically changing environment. We use the Avida Digital Evolution Platform (Ofria et al., 2009). Avida is an open-source system that has been used to conduct a wide range of well-regarded studies on evolutionary dynamics, including the origins of complex features (Lenski et al., 2003), the survival of the flattest effect (Wilke et al., 2001), and the origins of reproductive division of labor (Goldsby et al., 2014). Our experiments build directly on previous studies in Avida that characterized the de novo evolution of adaptive phenotypic plasticity (Clune et al., 2007; Lalejini and Ofria, 2016) as well as previous work investigating the evolutionary consequences of fluctuating environments for populations of non-plastic digital organisms (Li and Wilke, 2004; Canino-Koning et al., 2019). Of particular relevance, Clune et al. (2007) and Lalejini and Ofria (2016) experimentally demonstrated that adaptive phenotypic plasticity can evolve given the following four conditions (as described by Ghalambor et al., 2010): (1) populations experience temporal environmental variation, (2) these environments are differentiable by reliable cues, (3) each environment favors different phenotypic traits, and (4) no single phenotype exhibits high fitness across all environments. We build on this previous work, but we shift our focus from the evolutionary causes of adaptive phenotypic plasticity to investigate its evolutionary consequences in a fluctuating environment. Specifically, we examine the effects of adaptive plasticity on subsequent genomic and phenotypic change, the capacity to evolve and then maintain novel traits, and the accumulation of deleterious alleles.

Each of our experiments are divided into two phases: in phase one, we precondition sets of founder organisms with differing plastic or non-plastic adaptations; in phase two, we examine the subsequent evolution of populations founded with organisms from phase one under specific environmental conditions (Figure 2). First, we examine the evolutionary histories of phase two populations to test whether adaptive plasticity constrained subsequent genomic and phenotypic changes. Next, we evaluate how adaptive plasticity influences how well populations produced by each type of founder can evolve and retain novel adaptive traits. Finally, we examine lineages to determine whether adaptive plasticity facilitated the accumulation of cryptic genetic variation that would prove deleterious when the environment changed.

FIGURE 2

Figure 2. Overview of experimental design. The three plots in (A) show the environmental conditions used in every experiment and whether they reward or punish each base function. The two plots in (B) show the additional functions added in phases 2B and 2C. (C) Shows treatment differences and experimental phases. Treatments are listed on the left, with each treatment specifying its environmental configuration and whether sensors are functional. We conducted three independent two-phase experiments, each described on the right. Phases 2B and 2C are textured to match their function definitions in (B). Phase one is repeated for each experiment with 100 replicate populations per treatment per experiment. For each replicate at the end of phase one, we used an organism of the most abundant genotype to found the second phase population. All STATIC and NON-PLASTIC populations move on to phase two, but PLASTIC populations only continue to the second phase if their most abundant genotype exhibits optimal plasticity. Metrics are recorded only in phase two.

We found that the evolution of adaptive plasticity reduced subsequent rates of evolutionary change in a cyclic environment. The non-plastic populations underwent more frequent selective sweeps and accumulated many more genetic changes over time, as non-plastic populations relied on genetic variation from de novo mutations to continuously readapt to environmental changes. The evolution of adaptive phenotypic plasticity buffered populations against environmental fluctuations, whereas repeated selective sweeps in non-plastic populations drove the accumulation of deleterious mutations and the loss of secondary beneficial traits. As such, adaptively plastic populations were better able to retain novel traits than their non-plastic counterparts. In general, the evolution of adaptive phenotypic plasticity shifted evolutionary dynamics to be more similar to that of populations evolving in a static environment than to non-plastic populations evolving in an identical fluctuating environment.

2. Materials and Methods

2.1. The Avida Digital Evolution Platform

Avida is a study system wherein self-replicating computer programs (digital organisms) compete for space on a finite toroidal grid (Ofria et al., 2009). Each digital organism is defined by a linear sequence of program instructions (its genome) and a set of virtual hardware components used to interpret and express those instructions. Genomes are expressed sequentially except when the execution of one instruction (e.g., a “jump” instruction) deterministically changes which instruction should be executed next. Genomes are built using an instruction set that is both robust (i.e., any ordering of instructions is syntactically valid, though not necessarily meaningful) and Turing Complete (i.e., able to represent any computable function, though not necessarily in an efficient manner). The instruction set includes operations for basic computations, flow control (e.g., conditional logic and looping), input, output, and self-replication.

Organisms in Avida reproduce asexually by copying their genome instruction-by-instruction and then dividing. However, copy operations are imperfect and can result in single-instruction substitution mutations in an offspring's genome. For this work, we configured copy operations to err at a rate of one expected mutation for every 400 instructions copied (i.e., a per-instruction error rate of 0.0025). We held individual genomes at a fixed length of 100 instructions (the default genome size in Avida); that is, we did not include insertion and deletion mutations. We used fixed-length genomes to control for treatment-specific conditions resulting in the evolution of substantially different genome sizes (Lalejini and Ferguson, 2021a)¹, which could, on its own, drive differences in evolutionary outcomes among experimental treatments. When an organism divides in Avida, its offspring is placed in a random location on the toroidal grid, replacing any previous occupant. For this work, we used the default 60 × 60 grid size, which limits the maximum population size to 3,600 organisms. As such, improvements to the speed of self-replication are advantageous in the competition for space.

During evolution, organism replication rates improve in two ways: by improving genome efficiency (e.g., using a more compact encoding) or by accelerating the rate at which the genome is expressed (their “metabolic rate”). An organism's metabolic rate determines the speed at which it executes instructions in its genome. Initially, an organism's metabolic rate is proportional to the length of its genome, but that rate is adjusted as it completes designated functions, such as performing Boolean logic functions (Ofria et al., 2009). In this way, we can reward or punish particular phenotypic traits (i.e., Boolean logic functions).

2.1.1. Phenotypic Plasticity in Avida

In this work, we measure a digital organism's phenotype as the set of Boolean logic functions that it performs (i.e., expresses) in a given environment. Sensory instructions in the Avida instruction set allow organisms to detect how performing a particular function would affect their metabolic rate (see Supplementary Material for more details, Lalejini and Ferguson 2021a). We define a phenotypically plastic organism as one that uses sensory information to alter which functions it performs based on the environment.

Phenotypic plasticity in Avida can be adaptive or non-adaptive for a given set of environments. Adaptive plasticity shifts net function expression closer to the optimum for the given environments. Non-adaptive plasticity changes function expression in either a neutral or deleterious way. In this work, optimal plasticity toggles functions to always perfectly match the set of rewarded functions for the given set of environments.

2.2. Experimental Design

We conducted three independent experiments using Avida to investigate how the evolution of adaptive plasticity influences evolutionary outcomes in fluctuating environments. For each experiment, we compared the evolutionary outcomes of populations evolved under three treatments (Figure 2): (1) a PLASTIC treatment where the environment fluctuates, and digital organisms can use sensory instructions to differentiate between environmental states; (2) a NON-PLASTIC treatment with identical environment fluctuations, but where sensory instructions are disabled; and (3) a STATIC control where organisms evolve in a constant environment.

Each experiment was divided into two phases that each lasted for 200,000 updates² of evolution (Figure 2), which is equivalent to ~30,000–40,000 generations. In phase one of each experiment, we evolved plastic, non-plastic, and control organisms for use in phase two. In phase two, we founded new populations with these evolved organisms and examined their subsequent evolution under given combinations of treatment and experimental conditions. During phase two, we tracked and saved each population's evolutionary history as well as saving the full final population. Phase one was for preconditioning only; all comparisons between treatments were performed on phase two data.

2.2.1. Environments

We constructed three experimental environmental conditions, abbreviated hereafter as “ENV-A,” “ENV-B,” and “ENV-ALL.” In ENV-A, organisms are rewarded for expressing the NOT, AND, and OR Boolean logic functions, and organisms are punished for expressing the NAND, AND-NOT, and OR-NOT functions. ENV-B is the reverse of ENV-A; that is, in ENV-B, organisms are rewarded for expressing the NAND, AND-NOT, and OR-NOT functions and are punished for expressing the NOT, AND, and OR functions. In ENV-ALL, organisms are rewarded for expressing each of the NOT, AND, OR, NAND, AND-NOT, and OR-NOT functions. In all environmental conditions (ENV-A, ENV-B, and ENV-ALL), a rewarded function expressed by an organism doubles their metabolic rate, allowing them to execute twice as many instructions in the same amount of time. A punished function halves an organism's metabolic rate. Each Boolean logic function is a non-trivial trait to evolve, as they each require the coordination of multiple genetic instructions to express (Lenski et al., 2003).

In both the PLASTIC and NON-PLASTIC treatments, the environment cycles between equal-length periods of ENV-A and ENV-B. Each of these periods persist for 100 updates (~15–20 generations). Thus, populations experience a total of 1,000 full periods of ENV-A interlaced with 1,000 full periods of ENV-B during each experimental phase. Previous work has shown that this rate of change reliably allows for the evolution of adaptive phenotypic plasticity in Avida (Clune et al., 2007; Lalejini and Ofria, 2016).

Organisms in the PLASTIC treatments differentiate between ENV-A and ENV-B by executing one of six sensory instructions, each associated with a particular Boolean logic function; these sensory instructions detect whether their associated function is currently rewarded or punished. By using sensory information in combination with execution flow-control instructions, organisms can conditionally perform different functions depending on the current environmental conditions.

2.2.2. Experiment Phase 1—Environment Preconditioning

For each treatment, we founded 100 independent populations from a common ancestral strain capable only of self-replication. At the end of phase one, we identified the most abundant (i.e., dominant) genotype and sampled an organism with that genotype from each replicate population to found a new population for phase two.

For the PLASTIC treatment, we needed to ensure that our observations during the second phase of each experiment reflected the evolutionary consequences of adaptive plasticity. To do so, measured the plasticity of each PLASTIC-treatment population's dominant genotype by independently testing that genotype in each of ENV-A and ENV-B and recording the phenotype expressed in each environment. We discarded PLASTIC-treatment phase one populations if the dominant genotype did not exhibit optimal plasticity (as defined in section 2.1.1), which ensured that PLASTIC-treatment phase two populations were founded with an optimally plastic organism.

2.2.3. Experiment Phase 2A—Evolutionary Change Rate

We conducted experiment phase 2A to test for differences in evolutionary change—the accumulation of genetic and phenotypic changes—among populations evolving under each of our three treatment conditions (PLASTIC, NON-PLASTIC, and STATIC). Phase 2A continued exactly as phase one, except we tracked the rates of evolutionary change in each of the PLASTIC-, NON-PLASTIC-, and STATIC-treatment populations. Specifically, we quantified evolutionary change using four metrics (each described in Table 1): (1) coalescence event count, (2) mutation count, (3) phenotypic volatility, and (4) mutational robustness.

TABLE 1

Table 1. Metric descriptions.

While environmental conditions during phases one and 2A are identical, these phases are distinct in how populations are founded: each phase one population is founded with a common ancestor capable only of self-replication, whereas each phase two population is founded with an organism that was evolved in its treatment-specific conditions. Thus, during phase two, phylogenies are rooted with an ancestor that is well-adapted to its treatment conditions, which, in turn, ensures that our observations can exclude dynamics associated with initial adaptation.

2.2.4. Experiment Phase 2B—Novel Function Evolution

We conducted experiment phase 2B to quantify the extent to which organisms evolved under PLASTIC-, NON-PLASTIC-, and STATIC-treatment conditions were able to acquire and retain novel functions. Phase 2B extended the conditions of phase one by adding 71 novel Boolean logic functions (Ofria et al., 2009), which were always rewarded in all treatments. The original six phase one functions (NOT, NAND, AND, OR-NOT, OR, and AND-NOT; hereafter called “base” functions) continued to be rewarded or punished according to the particular treatment conditions. An organism's metabolic rate was increased by 10% for each novel function that it expressed (limited to one reward per function). This reward provided a selective pressure to evolve these functions, but their benefits did not overwhelm the 100% metabolic rate increase conferred by rewarded base functions. As such, populations in the PLASTIC and NON-PLASTIC treatments could not easily escape environmental fluctuations by abandoning the fluctuating base functions.

During this experiment, we used three metrics to quantify novel function acquisition and retention in evolving populations (each described in Table 1): (1) final novel function count, (2) novel function discovery, and (3) novel function loss.

2.2.5. Experiment Phase 2C—Deleterious Instruction Accumulation

We conducted experiment phase 2C to quantify the extent to which organisms evolved under PLASTIC-, NON-PLASTIC-, and STATIC-treatment conditions acquired deleterious instructions via mutation. Phase 2C extended the instruction set of phase one with a deleterious instruction. When an organism executes a deleterious instruction, it performs a “deleterious” function, which reduces the organism's metabolic rate (and thus reproductive success) but does not otherwise alter the organism's behavior. We imposed a 10% penalty each time an organism performed the deleterious function, making the deleterious instruction explicitly harmful to execute. We did not limit the number of times that an organism could perform the deleterious function, and as such, organisms could perform the deleterious function as many times as they executed the deleterious instruction.

We tracked the number of times each organism along the dominant lineage performed the deleterious function. Specifically, we used two metrics (each described in Table 1): (1) final deleterious function count and (2) deleterious function acquisition count.

2.3. Experimental Analyses

For each of our experiments, we tracked and analyzed the phylogenetic histories of evolving populations during phase two and generated a set of summary statistics (Table 1). For each phylogenetic history, we counted the number of times that the most recent common ancestor for the population shifted and used this value as the number of coalescence events. Next, at the final time point, we identified the most abundant genotype in the evolved population, and chose a representative organism from that genotype for further analysis. We used the lineage from the founding organism to the representative organism to summarize the evolutionary pathway of a population. This lineage represents the majority of the evolutionary history from the given population as long as the entire population traces back to the lineage in recent history. We manually inspected evolved phylogenies and found no evidence that any of our experimental treatments supported long-term coexistence. As such, analyses of the representative lineages reflect the majority of evolutionary history for a given population.

Some of our metrics (Table 1) required us to measure genotype-by-environment interactions. Importantly, in the fluctuating environments, we needed to differentiate phenotypic changes that were caused by mutations from those that were caused by environmental changes. To accomplish this, we produced organisms with the given focal genotype, measured their phenotype in each environmental condition, and aggregated the resulting phenotypes to create a phenotypic profile. Although organisms with different genotypes may express the same set of functions across environmental conditions, their phenotypic profiles may not necessarily be the same. For example, an organism that expresses NOT in ENV-A and NAND in ENV-B has a distinct phenotypic profile from one that expresses NAND in ENV-A and NOT in ENV-B. Because phenotypic profiles encapsulate function expression across all relevant environmental conditions (ENV-A and ENV-B), a change in phenotypic profile from parent to offspring indicates a mutationally-induced phenotypic change.

While most analyses employed here are retrospective metrics applied to lineages, digital evolution allows precise manipulations on individual organisms and genomes. Mutational robustness uses this technique when looking at the possible mutations on a representative genotype. Genomes in Avida are linear sequences of instructions, and as such, possible mutations can be simulated by substituting other instructions at the desired site. Indeed, the mutational robustness of a genotype examines all one-step mutations (i.e., each mutation where exactly one instruction is substituted). This allows us to disentangle whether the frequency of mutationally-induced phenotypic changes observed along a lineage is a consequence of evolved genetic architectures vs. the result of population dynamics that actively skewed the distribution of organisms along the lineage; that is, are genomes organized such that mutations are more likely to induce phenotypic changes, or are organisms with phenotypes different from their parents more likely to be successful and thus appear along the representative lineage?

2.4. Statistical Analyses

Across all of our experiments, we differentiated between sample distributions using non-parametric statistical tests. For each major analysis, we first performed a Kruskal-Wallis test (Kruskal and Wallis, 1952) to determine if there were significant differences in results from the PLASTIC, NON-PLASTIC, and STATIC treatments (significance level α = 0.05). If so, we applied a Wilcoxon rank-sum test (Wilcoxon, 1992) to distinguish between pairs of treatments. We applied Bonferroni corrections for multiple comparisons (Rice, 1989) where appropriate.

2.5. Software Availability

We conducted our experiments using a modified version of the Avida software, which is open source and freely available on GitHub (Lalejini and Ferguson, 2021a). Modifications to Avida included an improved phylogeny tracking system that enabled us to track coalescence events and the addition of custom sensory instructions specific to our experiments. We used Python for data processing, and we conducted all statistical analyses using R version 4 (R Core Team, 2021). We used the tidyverse collection of R packages (Wickham et al., 2019) to wrangle data, and we used the following R packages for analysis, graphing, and visualization: ggplot2 (Wickham et al., 2020), cowplot (Wilke, 2020), Color Brewer (Harrower and Brewer, 2003; Neuwirth, 2014), rstatix (Kassambara, 2021), ggsignif (Ahlmann-Eltze and Patil, 2021), scales (Wickham and Seidel, 2020), Hmisc (Harrell et al., 2020), fmsb (Nakazawa, 2019), and boot (Canty and Ripley, 2019). We used R markdown (Allaire et al., 2020) and bookdown (Xie, 2020) to generate web-enabled Supplementary Material. All of the source code for our experiments and analyses, including configuration files and guides for replication, can be found in our Supplementary Material, which is hosted on GitHub (Lalejini and Ferguson, 2021a). Additionally, our experimental data is available on the Open Science Framework at https://osf.io/sav2c/ (Lalejini and Ferguson, 2021b).

3. Results

3.1. Adaptive Phenotypic Plasticity Slows Evolutionary Change in Fluctuating Environments

In experimental phase 2A, we tested whether adaptive phenotypic plasticity constrained or promoted subsequent evolutionary change in a fluctuating environment. First, we compared the total amount of evolutionary change in populations evolved under the PLASTIC, NON-PLASTIC, and STATIC treatments as measured by coalescence event count, mutation count, and phenotypic volatility (Figure 3). According to each of these metrics, NON-PLASTIC populations experienced a larger magnitude of evolutionary change than either PLASTIC or STATIC populations. We observed significantly higher coalescence event counts in NON-PLASTIC populations than in PLASTIC or STATIC populations (Figure 3A). NON-PLASTIC lineages had significantly higher mutation counts (Figure 3B) and phenotypic volatility than PLASTIC or STATIC lineages (Figure 3C).

FIGURE 3

Figure 3. Magnitude of evolutionary change. Raincloud plots (Allen et al., 2019) of (A) coalescence event count, (B) mutation count, and (C) phenotypic volatility. See Table 1 for descriptions of each metric. Each plot is annotated with statistically significant comparisons (Bonferroni-corrected pairwise Wilcoxon rank-sum tests). Note that adaptive phenotypic plasticity evolved in 42 of 100 replicates from the PLASTIC treatment during phase one of this experiment; we used this more limited group to found 42 phase-two PLASTIC replicates from which we report these PLASTIC data.

Changing environments have been shown to increase generational turnover (i.e., how rapidly generations elapse) in Avida populations (Canino-Koning et al., 2016), which could explain why we observe a larger magnitude of evolutionary change at the end of 200,000 updates of evolution in NON-PLASTIC populations. Indeed, we found that significantly more generations of evolution elapsed in NON-PLASTIC populations (mean of 41, 090±2, 702 std. dev.) than in PLASTIC (mean of 31, 016±2, 615 std. dev.) or STATIC (mean of 30, 002±3, 011 std. dev.) populations during phase 2A (corrected Wilcoxon rank-sum tests, p < 10⁻⁴).

To evaluate whether increased generational turnover explains the greater magnitude of evolutionary change in NON-PLASTIC populations, we examined the average number of generations between coalescence events and the realized mutational robustness of lineages (Table 1). A coalescence event indicates a selective sweep, which is a hallmark of adaptive evolutionary change. Realized mutational robustness measures the frequency that mutations cause phenotypic changes along a lineage. We expect that static conditions should favor fit lineages with high realized mutational robustness that no longer undergo rapid adaptive change and hence do not trigger frequent coalescence events. Under fluctuating conditions, however, lineages must be composed of plastic organisms if they are to maintain both high fitness and realized mutational robustness. Without plasticity, we expect fluctuating conditions to produce lineages with low realized mutational robustness and frequent coalescence events as populations must continually acquire and fix mutations to readapt to the environment.

On average, significantly fewer generations elapsed between coalescence events in NON-PLASTIC populations than in either PLASTIC or STATIC populations (Figure 4A). We also found that both STATIC and PLASTIC lineages exhibited higher realized mutational robustness relative to that of NON-PLASTIC lineages (Figure 4B); that is, mutations observed along NON-PLASTIC lineages more often caused phenotypic changes in offspring. Overall, our results indicate that NON-PLASTIC populations underwent more rapid (and thus a greater amount of) evolutionary change than either PLASTIC or STATIC populations.

FIGURE 4

Figure 4. Pace of evolutionary change. Raincloud plots of (A) average number of generations between coalescence events, and (B) realized mutational robustness (Table 1). Each plot is annotated with statistically significant comparisons (Bonferroni-corrected pairwise Wilcoxon rank-sum tests).

While both STATIC and PLASTIC lineages exhibited high realized mutational robustness, we found that STATIC lineages exhibited higher realized robustness than PLASTIC lineages (Figure 4B). Overall, there were rare instances of mutations that caused a change in phenotypic profile across all PLASTIC lineages. Of these mutations, we found that over 80% (83 out of 102) of changes to phenotypic profiles were cryptic. That is, the mutations affected traits that would not have been expressed in the environment that the organism was born into but would have been expressed had the environment changed.

Given that NON-PLASTIC lineages exhibited the lowest realized mutational robustness of our three experimental treatments, we sought to determine if this effect was driven by differences in evolved genetic architectures. Specifically, did the NON-PLASTIC genetic architectures evolve such that mutations were more likely to result in phenotypic change? Such a mutational bias would trade off descendant fitness in the same environment in exchange for a chance of increasing descendant fitness in alternate environments. This strategy would be an example of diversifying bet-hedging (i.e., reducing expected mean fitness to lower variance in fitness; (Childs et al., 2010). Alternatively, the lower realized mutational robustness in NON-PLASTIC lineages could be due to survivorship bias, as we measured realized mutational robustness as the fraction of mutations observed along successful lineages that caused a phenotypic change. That is, our measure of realized mutational robustness concentrates only on mutations that appear in the representative lineage (i.e., “survivors” of selection), ignoring mutations that did not. Thus, lower realized mutational robustness in NON-PLASTIC lineages could simply be due to phenotype-altering mutations being more frequently favored by selection.

We analyzed the mutational robustness of representative genotypes by calculating the fraction of single-instruction mutations that change the phenotypic profile. We found that mutations to representative genotypes on NON-PLASTIC lineages are less likely to result in a phenotypic change than mutations to comparable genotypes on either STATIC or PLASTIC lineages (Figure 5). These data provide evidence against NON-PLASTIC lineages engaging in a mutation-driven bet-hedging strategy, and instead, are consistent with the hypothesis that lower realized mutational robustness in the NON-PLASTIC treatment was due to survivorship bias.

FIGURE 5

Figure 5. Mutational robustness. Raincloud plot of mutational robustness of each representative genotype (Table 1). The plot is annotated with statistically significant comparisons (Bonferroni-corrected pairwise Wilcoxon rank-sum tests).

In general, adaptive plasticity stabilized PLASTIC-treatment populations against environmental fluctuations, and their evolutionary dynamics more closely resembled those of populations evolving in a static environment. We observed no significant difference in the number and frequency of coalescence events in PLASTIC and STATIC populations. We did, however, observe small, but statistically significant, differences in each of the following metrics: elapsed generations, mutation counts, phenotypic volatility, realized mutational robustness, and mutational robustness between PLASTIC and STATIC populations. We expect that these differences are a result of plastic organisms needing to simultaneously maintain both function and function-regulation machinery, resulting in more genetic components that can be broken by mutations; moreover, many of these components are under relaxed selection in periods between environmental changes. Overall, these differences were not substantial enough to play an obvious role in any of the dynamics we analyzed, but could be examined further in the future.

3.2. Adaptively Plastic Populations Retain More Novel Functions Than Non-plastic Populations in Fluctuating Environments

We have so far shown that adaptive plasticity constrains the rate of evolutionary change in fluctuating environments. However, it is unclear how this dynamic influences the evolution of novel functions. Based on their relative rates of evolutionary change, we might expect NON-PLASTIC-treatment populations to evolve more novel functions than PLASTIC-treatment populations. But, how much of the evolutionary change in NON-PLASTIC populations is useful for exploring novel regions of the fitness landscape vs. continually rediscovering the same regions?

To answer this question, we quantified the number of novel functions performed by a representative organism in the final population of each replicate. We found that both PLASTIC and STATIC populations had significantly higher final function counts than NON-PLASTIC populations at the end of the experiment (Figure 6A). The final novel function count in PLASTIC and STATIC lineages could be higher than that of the NON-PLASTIC lineages for several non-mutually exclusive reasons. One possibility is that PLASTIC and STATIC lineages could be exploring a larger area of the fitness landscape when compared to NON-PLASTIC lineages. Another possibility is that the propensity of the NON-PLASTIC lineages to maintain novel traits could be significantly lower than PLASTIC or STATIC lineages. When we looked at the total sum of novel functions discovered by each of the PLASTIC, STATIC, and NON-PLASTIC lineages, we found that NON-PLASTIC lineages generally explored a larger area of the fitness landscape (Figure 6B). Although the NON-PLASTIC lineages discovered more novel functions, those lineages also exhibited significantly higher novel function loss when compared to PLASTIC and STATIC lineages (Figure 6C).

FIGURE 6

Figure 6. Novel function evolution. Raincloud plots of (A) final novel function count, (B) novel function discovery, and (C) novel function loss. See Table 1 for descriptions of each metric. Each plot is annotated with statistically significant comparisons (Bonferroni-corrected pairwise Wilcoxon rank-sum tests). Note that adaptive phenotypic plasticity evolved in 42 of 100 replicates from the PLASTIC treatment during phase one of this experiment; we used this more limited group to seed the resulting 42 phase-two PLASTIC replicates.

A larger number of generations elapsed in NON-PLASTIC populations than in PLASTIC or STATIC populations during our experiment (Lalejini and Ferguson, 2021a). Are NON-PLASTIC lineages discovering and losing novel functions more frequently than PLASTIC or STATIC lineages, or are our observations a result of differences in generational turnover? To answer this question, we converted the metrics of novel function discovery and novel function loss to rates by dividing each metric by the number of elapsed generations along the associated representative lineages. We found no significant difference in the frequency of novel function discovery between NON-PLASTIC and STATIC lineages, and we found that PLASTIC lineages had a lower frequency of novel function discovery than STATIC lineages (Figure 7A). Therefore, we cannot reject the possibility that the larger magnitude of function discovery in NON-PLASTIC lineages was driven by a larger number of elapsed generations. NON-PLASTIC lineages had a higher frequency of function loss than either PLASTIC or STATIC lineages, and PLASTIC lineages tended to have a lower frequency of novel function loss than STATIC lineages (Figure 7B).

FIGURE 7

Figure 7. Rate of novel function evolution. Raincloud plots of (A) novel function discovery frequency and (B) novel function loss frequency. Each plot is annotated with statistically significant comparisons (Bonferroni-corrected pairwise Wilcoxon rank-sum tests).

Next, we examined the frequency at which novel function loss along lineages co-occurred with the loss or gain of any of the six base functions. Across all NON-PLASTIC representative lineages, over 97% (10,998 out of 11,229) of instances of novel function loss co-occurred with a simultaneous change in base function profile. In contrast, across all PLASTIC and STATIC dominant lineages, we observed that approximately 20% (29 out of 142) and 2% (13 out of 631), respectively, of instances of novel function loss co-occurred with a simultaneous change in base function profile. As such, the losses of novel functions in NON-PLASTIC lineages appear to be primarily due to hitchhiking or epistatic effects where a mutation that knocks out a maladaptive base function (after the environment changes) also knocks out a beneficial novel function.

3.3. Lineages Without Plasticity That Evolve in Fluctuating Environments Express More Deleterious Functions

Phenotypic plasticity allows for genetic variation to accumulate in genomic regions that are unexpressed, which could lead to the fixation of deleterious instructions in PLASTIC populations. However, in NON-PLASTIC lineages, we observe a higher rate of novel function loss, indicating that they may be more susceptible to deleterious mutations (Figure 7B).

Therefore, in experiment phase 2C, we tested whether adaptive phenotypic plasticity can increase the incidence of deleterious function performance. Specifically, we added an instruction that triggered an explicitly deleterious function and measured the number of times it was executed. Each execution of the deleterious instruction reduces an organism's fitness by 10%. At the beginning of phase 2C, the deleterious instruction is not present in the population, as it was not part of the instruction set during phase one of evolution. Accordingly, if a deleterious instruction fixes in a population, it must be the result of evolutionary dynamics during phase 2C, including cryptic variation, hitchhiking, or as a result of sign epistasis where a deleterious instruction knocks out an even more maladaptive trait.

At the end of our experiment, no representative organisms from the PLASTIC or STATIC treatments performed the deleterious function under any environmental condition; however, representative organisms in 14% of replicates of the NON-PLASTIC treatment performed the deleterious function at least once. NON-PLASTIC lineages contained significantly more mutations that conferred the deleterious function as compared to PLASTIC or STATIC lineages (Figure 8A), and these mutations occurred at a significantly higher frequency in NON-PLASTIC lineages (Figure 8B).

FIGURE 8

Figure 8. Deleterious instruction accumulation. Raincloud plots of (A) deleterious function acquisition, (B) deleterious function acquisition frequency, and (C) the proportion of mutations that increase deleterious function expression along a lineage that co-occur with a change in phenotypic profile. Each plot is annotated with statistically significant comparisons (Bonferroni-corrected pairwise Wilcoxon rank-sum tests). Note that adaptive phenotypic plasticity evolved in 43 of 100 replicates from the PLASTIC treatment during phase one of this experiment; we used this more limited group to seed the 43 phase-two PLASTIC replicates.

Next, we measured how often mutations that increased deleterious function performance co-occurred with changes to the base function profile within representative lineages. A deleterious instruction can fix in a population by having a beneficial effect that outweighs its inherent cost (e.g., knocking out a punished function) or through linkage with a secondary beneficial mutation at another site within the genome. Across all NON-PLASTIC representative lineages, we found that ~49% (956 out of 1,916) of mutations that increased deleterious function expression co-occurred with a change in the base function profile (Figure 8C). In all representative lineages from the PLASTIC treatment, only 18 mutations increased deleterious function expression, and none co-occurred with a change in base function profile (Figure 8C). Likewise, only 58 mutations increased deleterious function performance in all representative lineages from the STATIC treatment, and none co-occurred with a change in base function profile (Figure 8C). We did not find compelling evidence that the few mutations that increased deleterious function expression occurred as cryptic variation in PLASTIC lineages.

We repeated this experiment with 3 and 30% metabolic rate penalties associated with the deleterious function, which produced results that were consistent with those reported here (Lalejini and Ferguson, 2021a).

4. Discussion

In this work, we used evolving populations of digital organisms to determine how adaptive phenotypic plasticity alters subsequent evolutionary dynamics and influences evolutionary outcomes in fluctuating environments. Specifically, we compared lineages of adaptively plastic organisms in fluctuating environments to both non-plastic organisms in those same environments and other non-plastic organisms in static environments.

4.1. Evolutionary Change

We found strong evidence that adaptive plasticity slows evolutionary change in fluctuating environments. Adaptively plastic populations experienced fewer coalescence events and fewer total genetic changes relative to non-plastic populations evolving under identical environmental conditions (Figure 3). Whereas non-plastic populations relied on de novo mutations to adapt to each environmental fluctuation, plastic populations leveraged sensory instructions to regulate function performance. Indeed, in fluctuating environments, selection pressures toggle after each environmental change. We hypothesize that in non-plastic populations such toggling would repeatedly drive the fixation of mutations that align an organism's phenotypic profile to the new conditions. This hypothesis is supported by the increased frequency of coalescence events in these populations (Figure 4A) as well as increased rates of genetic and phenotypic changes observed along the lineages of non-plastic organisms.

Representative lineages in the non-plastic treatment experienced lower realized mutational robustness than plastic and static lineages (Figure 4B). We reasoned that this lower realized mutational robustness was due to non-plastic populations evolving a bet-hedging strategy where mutations are more likely to modify the phenotypic profile. However, when we switched from measuring the realized mutational robustness of representative lineages to measuring the mutational robustness of representative genotypes (i.e., what fraction of one-step mutants change the phenotypic profile), we observed that non-plastic genotypes exhibited the highest mutational robustness of all three treatments (Figure 5). This result runs contrary to both our expectations and the results of other fluctuating environment studies in Avida (Canino-Koning et al., 2019). Canino-Koning et al. (2019) found that mutational robustness is negatively correlated with the number of function-encoding sites in the genome. In our work, most plastic and static genotypes encode all six base functions, while most non-plastic genotypes only encode functions from one environment; this results in fewer function-encoding sites, which may increase mutational robustness in non-plastic genotypes (relative to plastic and static genotypes). Regardless of the cause, this higher mutational robustness in non-plastic organisms indicates that bet-hedging is not driving the low realized mutational robustness observed in non-plastic lineages. Thus, we expect the lower realized mutational robustness in non-plastic lineages to be driven by survivorship bias. Because non-plastic lineages must rely on mutations to adapt to environmental changes, phenotype-altering mutations are often highly advantageous, and their selection decreases the realized mutational robustness of successful lineages.

To our knowledge, this study is the first in-depth empirical investigation into how the de novo evolution of adaptive plasticity shifts the course of subsequent evolution in a cyclic environment. The relative rates of evolutionary change that we observed in non-plastic populations, however, are consistent with results from previous digital evolution studies. For example, Dolson et al. (2020) showed that non-plastic populations that were evolved in cyclically changing environments exhibited higher phenotypic volatility and accumulated more mutations than that of populations evolved under static conditions. Furthermore, Lalejini and Ofria (2016) visually inspected the evolutionary histories of non-plastic organisms evolved in fluctuating environments, observing that mutations along successful lineages readily switched the set of traits expressed by offspring.

Our results are also consistent with conventional evolutionary theory. A trait's evolutionary response to selection depends on the strength of directional selection and on the amount of genetic variation for selection to act upon (Lande and Arnold, 1983; Zimmer and Emlen, 2013). In our experiments, non-plastic populations repeatedly experienced strong directional selection to toggle which functions were expressed after each environmental change. As such, retrospective analyses of successful lineages revealed rapid evolutionary responses (that is, high rates of genetic and phenotypic changes). Evolved adaptive plasticity shielded populations from strong directional selection when the environment changed by eliminating the need for a rapid evolutionary response to toggle function expression. Indeed, both theoretical and empirical studies have shown that adaptive plasticity can constrain evolutionary change by weakening directional selection on evolving populations (Price et al., 2003; Paenke et al., 2007; Ghalambor et al., 2015).

4.2. The Evolution and Maintenance of Novel Functions

In fluctuating environments, non-plastic populations explored a larger area of the fitness landscape than adaptively plastic populations (Figure 6B). However, adaptively plastic populations better exploited the fitness landscape, retaining a greater number of novel functions than non-plastic populations evolving under identical environmental conditions (Figure 6A). In our experiment, novel functions were less important to survival than the fluctuating base functions. In non-plastic populations, when a mutation changes a base function to better align with current environmental conditions, its benefit will often outweigh the cost of losing one or more novel functions. Indeed, we found that along non-plastic representative lineages, 97% of the mutations associated with novel function loss co-occurred with phenotypic changes that helped offspring adapt to current environmental conditions.

Previous studies have shown that transitory environmental changes can improve overall fitness landscape exploration in evolving populations of non-plastic digital organisms (Nahum et al., 2017). Similarly, changing environments have been shown to increase the rate of evolutionary adaptation in simulated network models (Kashtan et al., 2007). In our system, however, we found that repeated fluctuations reduced the ability of non-plastic populations to maintain and exploit functions; that said, we did find that repeated fluctuations may improve overall function discovery by increasing generational turnover. Consistent with our findings, Canino-Koning et al. (2019) found that non-plastic populations of digital organisms evolving in a cyclic environment maintained fewer novel traits than populations evolving in static environments.

Our results suggest that adaptive phenotypic plasticity can improve the potential for populations to exploit novel resources by stabilizing them against stressful environmental changes. The stability that we observe may also lend some support to the hypothesis that phenotypic plasticity can rescue populations from extinction under changing environmental conditions (Chevin et al., 2010).

Our data do not necessarily provide evidence for or against the genes as followers hypothesis. The genes as followers hypothesis focuses on contexts where plastic populations experience novel or abnormally stressful environmental change. However, in our system, environmental changes were cyclic (not novel), and no single environmental change was abnormally stressful. Further, the introduction of novel functions during the second phase of the experiment merely added static opportunities for fitness improvement. This addition did not change the meaning of existing environmental cues, nor did it require those cues to be used in new ways.

4.3. The Accumulation of Deleterious Alleles

We found that non-plastic lineages that evolved in a fluctuating environment exhibited both greater totals and higher rates of deleterious function acquisition than that of adaptively plastic lineages (Figure 8). There are several, non-mutually exclusive possibilities that could explain the fixation of explicitly deleterious instructions: random genetic drift, deleterious hitchhiking, epistatic effects, and cryptic variation (in plastic organisms). We find it unlikely that random genetic drift explains our observations. Each time an organism expresses a deleterious instruction, the organism incurs a 10% penalty to their replication rate, which results in strong purifying selection against mutations that cause offspring to execute deleterious instructions.

In asexual populations without horizontal gene transfer, all co-occurring mutations are linked. As such, deleterious mutations linked with a stronger beneficial mutation (i.e., a driver) can sometimes “hitchhike” to fixation (Smith and Haigh, 1974; Buskirk et al., 2017; Van den Bergh et al., 2018). Natural selection normally prevents deleterious mutations from reaching high frequencies, as such mutants are outcompeted. However, when a beneficial mutation sweeps to fixation in a clonal population, it carries along any linked genetic material, including other beneficial, neutral, or deleterious mutations (Smith and Haigh, 1974; Barton, 2000). Therefore, deleterious genetic hitchhiking could have contributed to deleterious instruction accumulation along non-plastic lineages in changing environments.

Epistatic effects (i.e., interactions between genes) could have also contributed to deleterious instruction accumulation along non-plastic lineages. On their own, mutations that increase deleterious instruction execution are maladaptive; however, if such a mutation were to also knock out an even more harmful function, that mutation may have a net beneficial effect. As such, mutations that confer increased deleterious instruction execution could directly drive a selective sweep.

Representative lineages from non-plastic populations in the cyclic environment exhibited higher mutation accumulation (Figure 3B), novel function loss (Figure 6C), and deleterious function acquisition (Figure 8A) than their plastic counterparts. In aggregate, we found that many (~49%; 956/1,916) mutations that increased deleterious instruction execution in offspring co-occurred with mutations that provided an even stronger benefit by adapting the offspring to an environmental change. We expect that an even larger fraction of these deleterious mutations were linked to beneficial mutations, but our analysis only counted mutations that co-occurred in the same generation. Our analyses did not distinguish between epistatic effects and deleterious hitchhiking; however, more fine-grained analyses of secondary effects of mutations that conferred deleterious instruction execution could be performed in future work to disentangle these two mechanisms.

Theory predicts that under relaxed selection deleterious mutations should accumulate as cryptic variation in unexpressed traits (Lahti et al., 2009). Contrary to this expectation, we did not find evidence of deleterious instructions accumulating as cryptic variation in adaptively plastic lineages. One possible explanation is that the period of time between environmental changes was too brief for variants carrying unexpressed deleterious instructions to drift to high frequencies before the environment changed, after which purifying selection would have removed such variants. Indeed, we would not expect drift to fix an unexpressed trait since we tuned the frequency of environmental fluctuations to prevent valuable traits from being randomly eliminated during the off environment. Additionally, plastic organisms in Avida usually adjust their phenotype by toggling the expression of a minimal number of key instructions, leaving little genomic space for cryptic variation to accumulate.

4.4. Limitations and Future Directions

Our work lays the groundwork for using digital evolution experiments to investigate the evolutionary consequences of phenotypic plasticity in a range of contexts. However, the data presented here are limited to the evolution of adaptively plastic populations. Future work might explore the evolutionary consequences of maladaptive and non-adaptive phenotypic plasticity (e.g., Leroi et al. 1994), which are known to bias evolutionary outcomes (Ghalambor et al., 2015).

Additionally in our experiments, sensory instructions perfectly differentiated between ENV-A and ENV-B, and environmental fluctuations never exposed populations to entirely new conditions. These parameters have been shown to influence evolutionary outcomes (Li and Wilke, 2004; Boyer et al., 2021), which if relaxed in the context of further digital evolution experiments, may yield additional insights. Our experiments also focused on asexually reproducing digital organisms. Sexual reproduction has been shown to be advantageous in rapidly changing environments, such as the cyclic environments used in our study (Misevic et al., 2010). Future work could investigate how sexual reproduction affects the evolutionary consequences of adaptive plasticity.

We focused our analyses on the lineages of organisms with the most abundant genotype in the final population. These successful lineages represented the majority of the evolutionary histories of populations at the end of our experiment, as populations did not exhibit long-term coexistence of different clades. Our analyses, therefore, gave us an accurate picture of what fixed in the population. We did not, however, examine the lineages of extinct clades. Future work will extend our analyses to include extinct lineages, giving us a more complete view of evolutionary history, which may allow us to better distinguish adaptively plastic populations from populations evolving in a static environment.

As with any wet-lab experiment, our results are in the context of a particular model organism: “Avidian” self-replicating computer programs. Digital organisms in Avida regulate responses to environmental cues using a combination of sensory instructions and conditional logic instructions (if statements). The if instructions conditionally execute a single instruction depending on previous computations and the state of memory. As such, plastic organisms in Avida typically regulate phenotypes by toggling the expression of a small number of key instructions as opposed to regulating cohorts of instructions under the control of a single regulatory sequence (Lalejini and Ferguson, 2021a). This bias may limit the accumulation of hidden genetic variation in Avida genomes. However, as there are many model biological organisms, there are many model digital organisms that have different regulatory mechanisms (e.g., Lalejini and Ofria 2018) that should be used to test the generality of our results.

As with most digital evolution experiments, our mutation rates were high and population sizes were small (3,600 individuals) relative to experiments with microbes or conditions common in nature. As such, beneficial mutations can be generated rapidly and selective sweeps can occur quickly. Moreover, our analyses were limited to a single rate of environmental change and simple function reward structures, which likely influenced the rates of selective sweeps observed in our experiments. Future studies could address these limitations by increasing population sizes, decreasing mutation rates, investigating different function rewards and punishments, and altering the time spent in each environment.

Data Availability Statement

The datasets generated and analyzed for this study can be found on the Open Science Framework at https://osf.io/sav2c/ (Lalejini and Ferguson, 2021b).

Author Contributions

AL and AF designed the experiments, developed the necessary experiment software, conducted experiments, analyzed the results, and drafted the manuscript. AL, AF, NG, and CO edited and approved the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Department of Energy through Grant No. DE-SC0019436 and by the National Science Foundation (NSF) through the BEACON Center (DBI-0939454), a Graduate Research Fellowship to AL (DGE-1424871), and NSF Grant No. DEB-1655715.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank members of the MSU Digital Evolution Lab for helpful comments and suggestions on this work. We also thank Luis Zaman whose keynote talk at the 2020 Artificial Life conference inspired the initial exploratory experiments from which this work blossomed. Additionally, we thank the reviewers whose thoughtful comments and suggestions improved the quality of our manuscript. MSU provided computational resources through the Institute for Cyber-Enabled Research. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation, Department of Energy, Michigan State University, or The University of Idaho.

Supplementary Material

The Supplementary Material for this article is hosted on GitHub and can be found online at https://github.com/amlalejini/evolutionary-consequences-of-plasticity (Lalejini and Ferguson, 2021a).

Footnotes

1. ^We repeated our experiments without genome size restrictions and observed qualitatively similar results (see Supplementary Material, Lalejini and Ferguson 2021a).

2. ^One update in Avida is the amount of time required for the average organism to execute 30 instructions. See Ofria et al. (2009) for more details.

References

Ahlmann-Eltze, C., and Patil, I. (2021). ggsignif: Significance Brackets for 'ggplot2'. R Package Version 0.6.1.

Allaire, J., Xie, Y., McPherson, J., Luraschi, J., Ushey, K., Atkins, A., et al. (2020). rmarkdown: Dynamic Documents for R. R Package Version 2.6.

Allen, M., Poggiali, D., Whitaker, K., Marshall, T. R., and Kievit, R. A. (2019). Raincloud plots: a multi-platform tool for robust data visualization. Wellcome Open Res. 4:63. doi: 10.12688/wellcomeopenres.15191.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ancel, L. W. (2000). Undermining the Baldwin expediting effect: does phenotypic plasticity accelerate evolution? Theor. Popul. Biol. 58, 307–319. doi: 10.1006/tpbi.2000.1484

PubMed Abstract | CrossRef Full Text | Google Scholar

Barrett, R., and Schluter, D. (2008). Adaptation from standing genetic variation. Trends Ecol. Evol. 23, 38–44. doi: 10.1016/j.tree.2007.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Barton, N. H. (2000). Genetic hitchhiking. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 355, 1553–1562. doi: 10.1098/rstb.2000.0716

PubMed Abstract | CrossRef Full Text | Google Scholar

Beaumont, H. J. E., Gallie, J., Kost, C., Ferguson, G. C., and Rainey, P. B. (2009). Experimental evolution of bet hedging. Nature 462, 90–93. doi: 10.1038/nature08504

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyer, S., Hérissant, L., and Sherlock, G. (2021). Adaptation is influenced by the complexity of environmental change during evolution in a dynamic environment. PLoS Genet. 17:e1009314. doi: 10.1371/journal.pgen.1009314

PubMed Abstract | CrossRef Full Text | Google Scholar

Buskirk, S. W., Peace, R. E., and Lang, G. I. (2017). Hitchhiking and epistasis give rise to cohort dynamics in adapting populations. Proc. Natl. Acad. Sci. U.S.A. 114, 8330–8335. doi: 10.1073/pnas.1702314114

PubMed Abstract | CrossRef Full Text | Google Scholar

Canino-Koning, R., Wiser, M. J., and Ofria, C. (2016). “The evolution of evolvability: changing environments promote rapid adaptation in digital organisms,” in Proceedings of the Artificial Life Conference 2016 (Cancun: MIT Press), 268–275. doi: 10.7551/978-0-262-33936-0-ch047

CrossRef Full Text | Google Scholar

Canino-Koning, R., Wiser, M. J., and Ofria, C. (2019). Fluctuating environments select for short-term phenotypic variation leading to long-term exploration. PLoS Comput. Biol. 15:e1006445. doi: 10.1371/journal.pcbi.1006445

PubMed Abstract | CrossRef Full Text | Google Scholar

Canty, A., and Ripley, B. D. (2019). boot: Bootstrap R (S-Plus) Functions. R Package Version 1.3-23.

Chadwick, W., and Little, T. J. (2005). A parasite-mediated life-history shift in Daphnia magna. Proc. R. Soc. B Biol. Sci. 272, 505–509. doi: 10.1098/rspb.2004.2959

PubMed Abstract | CrossRef Full Text | Google Scholar

Chevin, L.-M., and Lande, R. (2010). When do adaptive plasticity and genetic evolution prevent extinction of a density-regulated population? Evolution 64, 1143–1150. doi: 10.1111/j.1558-5646.2009.00875.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chevin, L.-M., Lande, R., and Mace, G. M. (2010). Adaptation, plasticity, and extinction in a changing environment: towards a predictive theory. PLoS Biol. 8:e1000357. doi: 10.1371/journal.pbio.1000357

PubMed Abstract | CrossRef Full Text | Google Scholar

Childs, D. Z., Metcalf, C., and Rees, M. (2010). Evolutionary bet-hedging in the real world: empirical evidence and challenges revealed by plants. Proc. R. Soc. B Biol. Sci. 277, 3055–3064. doi: 10.1098/rspb.2010.0707

PubMed Abstract | CrossRef Full Text | Google Scholar

Clune, J., Ofria, C., and Pennock, R. T. (2007). “Investigating the emergence of phenotypic plasticity in evolving digital organisms,” in Advances in Artificial Life, Vol. 4648, eds F. Almeida e Costa, L. M. Rocha, E. Costa, I. Harvey, and A. Coutinho (Berlin; Heidelberg: Springer Berlin Heidelberg), 74–83. doi: 10.1007/978-3-540-74913-4_8

CrossRef Full Text | Google Scholar

Covert, A. W., Lenski, R. E., Wilke, C. O., and Ofria, C. (2013). Experiments on the role of deleterious mutations as stepping stones in adaptive evolution. Proc. Natl. Acad. Sci. U.S.A. 110, E3171–E3178. doi: 10.1073/pnas.1313424110

PubMed Abstract | CrossRef Full Text | Google Scholar

Crispo, E. (2007). The Baldwin effect and genetic assimilation: revisiting two mechanisms of evolutionary change mediated by phenotypic plasticity. Evolution 61, 2469–2479. doi: 10.1111/j.1558-5646.2007.00203.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Dolson, E., Lalejini, A., Jorgensen, S., and Ofria, C. (2020). Interpreting the tape of life: ancestry-based analyses provide insights and intuition about evolutionary dynamics. Artif. Life 26, 58–79. doi: 10.1162/artl_a_00313

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunn, A. M., Hogg, J. C., Kelly, A., and Hatcher, M. J. (2005). Two cues for sex determination in Gammarus duebeni: adaptive variation in environmental sex determination? Limnol. Oceanogr. 50, 346–353. doi: 10.4319/lo.2005.50.1.0346

CrossRef Full Text | Google Scholar

Forsman, A. (2015). Rethinking phenotypic plasticity and its consequences for individuals, populations and species. Heredity 115, 276–284. doi: 10.1038/hdy.2014.92

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghalambor, C. K., Angeloni, L. M., and Carroll, S. P. (2010). “Behavior as phenotypic plasticity,” in Evolutionary Behavioral Ecology, eds D. Westneat and C. W. Fox (New York, NY: Oxford University Press), 90–107.

Google Scholar

Ghalambor, C. K., Hoke, K. L., Ruell, E. W., Fischer, E. K., Reznick, D. N., and Hughes, K. A. (2015). Non-adaptive plasticity potentiates rapid adaptive evolution of gene expression in nature. Nature 525, 372–375. doi: 10.1038/nature15256

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghalambor, C. K., McKay, J. K., Carroll, S. P., and Reznick, D. N. (2007). Adaptive versus non-adaptive phenotypic plasticity and the potential for contemporary adaptation in new environments. Funct. Ecol. 21, 394–407. doi: 10.1111/j.1365-2435.2007.01283.x

CrossRef Full Text | Google Scholar

Gibert, P., Debat, V., and Ghalambor, C. K. (2019). Phenotypic plasticity, global change, and the speed of adaptive evolution. Curr. Opin. Insect Sci. 35, 34–40. doi: 10.1016/j.cois.2019.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gibson, G., and Dworkin, I. (2004). Uncovering cryptic genetic variation. Nat. Rev. Genet. 5, 681–690. doi: 10.1038/nrg1426

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldsby, H. J., Knoester, D. B., Ofria, C., and Kerr, B. (2014). The evolutionary origin of somatic cells under the dirty work hypothesis. PLoS Biol. 12:e1001858. doi: 10.1371/journal.pbio.1001858

PubMed Abstract | CrossRef Full Text | Google Scholar

Gupta, A. P., and Lewontin, R. C. (1982). A study of reaction norms in natural populations of Drosophila pseudoobscura. Evolution 36:934. doi: 10.2307/2408073

CrossRef Full Text | Google Scholar

Harrell F. E. Jr, with contributions from Charles Dupont, and many others. (2020). HMISC: Harrell Miscellaneous. R Package Version 4.4-2.

Harrower, M., and Brewer, C. A. (2003). ColorBrewer.org: an online tool for selecting colour schemes for maps. Cartogr. J. 40, 27–37. doi: 10.1179/000870403235002042

CrossRef Full Text | Google Scholar

Hendry, A. P. (2016). Key questions on the role of phenotypic plasticity in eco-evolutionary dynamics. J. Hered. 107, 25–41. doi: 10.1093/jhered/esv060

PubMed Abstract | CrossRef Full Text | Google Scholar

Huey, R., Hertz, P., and Sinervo, B. (2003). Behavioral drive versus behavioral inertia in evolution: a null model approach. Am. Nat. 161, 357–366. doi: 10.1086/346135

PubMed Abstract | CrossRef Full Text | Google Scholar

Kashtan, N., Noor, E., and Alon, U. (2007). Varying environments can speed up evolution. Proc. Natl. Acad. Sci. U.S.A. 104, 13711–13716. doi: 10.1073/pnas.0611630104

PubMed Abstract | CrossRef Full Text | Google Scholar

Kassambara, A. (2021). rstatix: Pipe-Friendly Framework for Basic Statistical Tests. R Package Version 0.7.0.

Kruskal, W. H., and Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621. doi: 10.1080/01621459.1952.10483441

CrossRef Full Text | Google Scholar

Lahti, D. C., Johnson, N. A., Ajie, B. C., Otto, S. P., Hendry, A. P., Blumstein, D. T., et al. (2009). Relaxed selection in the wild. Trends Ecol. Evol. 24, 487–496. doi: 10.1016/j.tree.2009.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Lalejini, A., and Ferguson, A. (2021a). Supplemental Material. GitHub. doi: 10.5281/zenodo.4642704

CrossRef Full Text | Google Scholar

Lalejini, A., and Ferguson, A. J. (2021b). Data for Adaptive Phenotypic Plasticity Stabilizes Evolution in Fluctuating Environments. OSF. Available online at: https://osf.io/sav2c.

Google Scholar

Lalejini, A., and Ofria, C. (2016). “The evolutionary origins of phenotypic plasticity,” in Proceedings of the Artificial Life Conference 2016 (Cancun: MIT Press), 372–379. doi: 10.7551/978-0-262-33936-0-ch063

PubMed Abstract | CrossRef Full Text | Google Scholar

Lalejini, A., and Ofria, C. (2018). “Evolving reactive agents with SignalGP,” in The 2018 Conference on Artificial Life (Tokyo: MIT Press), 368–369. doi: 10.1162/isal_a_00069

CrossRef Full Text | Google Scholar

Lande, R., and Arnold, S. J. (1983). The measurement of selection on correlated characters. Evolution 37:1210. doi: 10.2307/2408842

CrossRef Full Text | Google Scholar

Lenski, R. E., Ofria, C., Collier, T. C., and Adami, C. (1999). Genome complexity, robustness and genetic interactions in digital organisms. Nature 400, 661–664. doi: 10.1038/23245

PubMed Abstract | CrossRef Full Text | Google Scholar

Lenski, R. E., Ofria, C., Pennock, R. T., and Adami, C. (2003). The evolutionary origin of complex features. Nature 423, 139–144. doi: 10.1038/nature01568

PubMed Abstract | CrossRef Full Text | Google Scholar

Leroi, A. M., Bennett, A. F., and Lenski, R. E. (1994). Temperature acclimation and competitive fitness: an experimental test of the beneficial acclimation assumption. Proc. Natl. Acad. Sci. U.S.A. 91, 1917–1921. doi: 10.1073/pnas.91.5.1917

PubMed Abstract | CrossRef Full Text | Google Scholar

Levis, N. A., and Pfennig, D. W. (2016). Evaluating “plasticity-first” evolution in nature: key criteria and empirical approaches. Trends Ecol. Evol. 31, 563–574. doi: 10.1016/j.tree.2016.03.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., and Wilke, C. O. (2004). Digital evolution in time-dependent fitness landscapes. Artif. Life 10, 123–134. doi: 10.1162/106454604773563559

PubMed Abstract | CrossRef Full Text | Google Scholar

Misevic, D., Ofria, C., and Lenski, R. E. (2010). Experiments with digital organisms on the origin and maintenance of sex in changing environments. J. Hered. 101(Suppl. 1), S46–S54. doi: 10.1093/jhered/esq017

PubMed Abstract | CrossRef Full Text | Google Scholar

Nahum, J. R., West, J., Althouse, B. M., Zaman, L., Ofria, C., and Kerr, B. (2017). “Improved adaptation in exogenously and endogenously changing environments,” in Proceedings of the 14th European Conference on Artificial Life ECAL 2017 (Lyon: MIT Press), 306–313. doi: 10.7551/ecal_a_052

CrossRef Full Text | Google Scholar

Nakazawa, M. (2019). fmsb: Functions for Medical Statistics Book With Some Demographic Data. R Package Version 0.7.0.

Neuwirth, E. (2014). RColorBrewer: ColorBrewer Palettes. R Package Version 1.1-2.

Ofria, C., Bryson, D. M., and Wilke, C. O. (2009). “Avida: a software platform for research in computational evolutionary biology,” in Artificial Life Models in Software, eds M. Komosinski and A. Adamatzky (London: Springer London), 3–35. doi: 10.1007/978-1-84882-285-6_1

PubMed Abstract | CrossRef Full Text | Google Scholar

Paaby, A. B., and Rockman, M. V. (2014). Cryptic genetic variation: evolution's hidden substrate. Nat. Rev. Genet. 15, 247–258. doi: 10.1038/nrg3688

PubMed Abstract | CrossRef Full Text | Google Scholar

Paenke, I., Sendhoff, B., and Kawecki, T. (2007). Influence of plasticity and learning on evolution under directional selection. Am. Nat. 170, E47–E58. doi: 10.1086/518952

PubMed Abstract | CrossRef Full Text | Google Scholar

Pigliucci, M. (2006). Phenotypic plasticity and evolution by genetic assimilation. J. Exp. Biol. 209, 2362–2367. doi: 10.1242/jeb.02070

PubMed Abstract | CrossRef Full Text | Google Scholar

Price, T. D., Qvarnstrm, A., and Irwin, D. E. (2003). The role of phenotypic plasticity in driving genetic evolution. Proc. R. Soc. Lond. Ser. B Biol. Sci. 270, 1433–1440. doi: 10.1098/rspb.2003.2372

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team (2021). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.

Google Scholar

Rice, W. R. (1989). Analyzing tables of statistical tests. Evolution 43:223. doi: 10.1111/j.1558-5646.1989.tb04220.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Schaum, C. E., and Collins, S. (2014). Plasticity predicts evolution in a marine alga. Proc. R. Soc. B Biol. Sci. 281:20141486. doi: 10.1098/rspb.2014.1486

PubMed Abstract | CrossRef Full Text | Google Scholar

Schlichting, C. D. (2008). Hidden reaction norms, cryptic genetic variation, and evolvability. Ann. N. Y. Acad. Sci. 1133, 187–203. doi: 10.1196/annals.1438.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Schlichting, C. D., and Wund, M. A. (2014). Phenotypic plasticity and epigenetic marking: an assessment of evidence for genetic accommodation. Evolution 68, 656–672. doi: 10.1111/evo.12348

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwander, T., and Leimar, O. (2011). Genes as leaders and followers in evolution. Trends Ecol. Evol. 26, 143–151. doi: 10.1016/j.tree.2010.12.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, J. M., and Haigh, J. (1974). The hitch-hiking effect of a favourable gene. Genet. Res. 23, 23–35. doi: 10.1017/S0016672300014634

CrossRef Full Text | Google Scholar

Van den Bergh, B., Swings, T., Fauvart, M., and Michiels, J. (2018). Experimental design, population dynamics, and diversity in microbial experimental evolution. Microbiol. Mol. Biol. Rev. 82:e00008–18. doi: 10.1128/MMBR.00008-18

PubMed Abstract | CrossRef Full Text | Google Scholar

Wennersten, L., and Forsman, A. (2012). Population-level consequences of polymorphism, plasticity and randomized phenotype switching: a review of predictions. Biol. Rev. 87, 756–767. doi: 10.1111/j.1469-185X.2012.00231.x

PubMed Abstract | CrossRef Full Text | Google Scholar

West-Eberhard, M. J. (2003). Developmental Plasticity and Evolution. (New York, NY: Oxford University Press). doi: 10.1093/oso/9780195122343.001.0001

CrossRef Full Text | Google Scholar

West-Eberhard, M. J. (2005). Developmental plasticity and the origin of species differences. Proc. Natl. Acad. Sci. U.S.A. 102(Suppl. 1), 6543–6549. doi: 10.1073/pnas.0501844102

PubMed Abstract | CrossRef Full Text | Google Scholar

West-Eberhard, M. J. (2008). “Phenotypic plasticity,” in Encyclopedia of Ecology, eds S. E. Jørgensen and B. D. Fath (Amsterdam: Elsevier), 2701–2707. doi: 10.1016/B978-008045405-4.00837-5

CrossRef Full Text | Google Scholar

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., Franois, R., et al. (2019). Welcome to the tidyverse. J. Open Source Softw. 4:1686. doi: 10.21105/joss.01686

CrossRef Full Text | Google Scholar

Wickham, H., Chang, W., Henry, L., Pedersen, T. L., Takahashi, K., Wilke, C., et al. (2020). ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. R Package Version 3.3.2. (New York, NY).

Wickham, H., and Seidel, D. (2020). Scales: Scale Functions for Visualization. R Package Version 1.1.1

Wilcoxon, F. (1992). “Individual comparisons by ranking methods,” in Breakthroughs in Statistics, eds S. Kotz and N. L. Johnson (New York, NY: Springer New York), 196–202. doi: 10.1007/978-1-4612-4380-9_16

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilke, C. O. (2020). cowplot: Streamlined Plot Theme and Plot Annotations for ggplot2. R Package Version 1.1.0.

Wilke, C. O., and Adami, C. (2002). The biology of digital organisms. Trends Ecol. Evol. 17, 528–532. doi: 10.1016/S0169-5347(02)02612-5

CrossRef Full Text | Google Scholar

Wilke, C. O., Wang, J. L., Ofria, C., Lenski, R. E., and Adami, C. (2001). Evolution of digital organisms at high mutation rates leads to survival of the flattest. Nature 412, 331–333. doi: 10.1038/35085569

PubMed Abstract | CrossRef Full Text | Google Scholar

Winger, B. M., Auteri, G. G., Pegan, T. M., and Weeks, B. C. (2019). A long winter for the Red Queen: rethinking the evolution of seasonal migration. Biol. Rev. 94, 737–752. doi: 10.1111/brv.12476

PubMed Abstract | CrossRef Full Text | Google Scholar

Wund, M. A. (2012). Assessing the impacts of phenotypic plasticity on evolution. Integr. Compar. Biol. 52, 5–15. doi: 10.1093/icb/ics050

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y. (2020). bookdown: Authoring Books and Technical Documents With R Markdown. R Package Version 0.21.

Google Scholar

Zheng, J., Payne, J. L., and Wagner, A. (2019). Cryptic genetic variation accelerates evolution by opening access to diverse adaptive peaks. Science 365, 347–353. doi: 10.1126/science.aax1837

PubMed Abstract | CrossRef Full Text | Google Scholar

Zimmer, C., and Emlen, D. J. (2013). Evolution: Making Sense of Life. Greenwood Village, CO: Roberts and Company Publishers.