
==== Front
G3 (Bethesda)
Genetics
g3journal
G3: Genes | Genomes | Genetics
2160-1836
Oxford University Press US

38984708
10.1093/g3journal/jkae144
jkae144
Investigation
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
On the evolutionary origin of discrete phenotypic plasticity
https://orcid.org/0000-0002-0830-2009
Sakamoto Takahiro SOKENDAI, Research Center for Integrative Evolutionary Science, The Graduate University for Advanced Studies, Hayama, Kanagawa 240-0193, Japan
Department of Genomics and Evolutionary Biology, National Institute of Genetics, 1111 Yata, Mishima, Shizuoka 411-8540, Japan

https://orcid.org/0000-0001-6375-9231
Innan Hideki SOKENDAI, Research Center for Integrative Evolutionary Science, The Graduate University for Advanced Studies, Hayama, Kanagawa 240-0193, Japan

Jakobsson M Editor
Corresponding author: SOKENDAI, Research Center for Integrative Evolutionary Science, The Graduate University for Advanced Studies, Hayama, Kanagawa 240-0193, Japan. Email: innan_hideki@soken.ac.jp
Corresponding author: SOKENDAI, Research Center for Integrative Evolutionary Science, The Graduate University for Advanced Studies, Hayama, Kanagawa 240-0193, Japan. Email: sakamoto_takahiro@nig.ac.jp
Conflicts of interest The authors declare no competing interest.

9 2024
10 7 2024
10 7 2024
14 9 jkae14403 5 2024
24 6 2024
19 7 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of The Genetics Society of America.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Phenotypic plasticity provides an attractive strategy for adapting to various environments, but the evolutionary mechanism of the underlying genetic system is poorly understood. We use a simple gene regulatory network model to explore how a species acquires phenotypic plasticity, particularly focusing on discrete phenotypic plasticity, which has been difficult to explain by quantitative genetic models. Our approach employs a population genetic framework that integrates the developmental process, where each individual undergoes growth to develop its phenotype, which subsequently becomes subject to selection pressures. Our model considers two alternative types of environments, with the gene regulatory network including a sensor gene that turns on and off depending on the type of environment. With this assumption, we demonstrate that the system gradually adapts by acquiring the ability to produce two distinct optimum phenotypes under two types of environments without changing genotype, resulting in phenotypic plasticity. We find that the resulting plasticity is often discrete after a lengthy period of evolution. Our results suggest that gene regulatory networks have a notable capacity to flexibly produce various phenotypes in response to environmental changes. This study also shows that the evolutionary dynamics of phenotype may differ significantly between mechanistic-based developmental models and quantitative genetics models, suggesting the utility of incorporating gene regulatory networks into evolutionary models.

phenotypic plasticity
population genetics
development
gene regulatory network
Graduate University for Advanced Studies, SOKENDAI Japan Society for the Promotion of Science 10.13039/501100001691 JSPS KAKENHI 23H02521 23H04736 JSPS KAKENHI JP20J21760 JP23KJ2158
==== Body
pmcIntroduction

Phenotypic plasticity is the ability to produce different phenotypes from the same genotype depending on the environment, which can be an attractive strategy in adaptation. Such plasticity should confer a tremendous selective advantage in temporally or spatially fluctuating environments when each phenotype is adaptive in its specific environment. Of special interest is the cases where multiple different phenotypes arise depending on the environment with almost no intermediate phenotypes between them (i.e. discrete phenotypic plasticity), that is, the distribution of phenotype is intrinsically discrete (Fig. 1a). This situation is quite different from continuous plasticity, in which phenotypic trait changes gradually as environment changes (Fig. 1b). There are a number of example of continuous plasticity observed in nature (Huey and Kingsolver 1989; Brakefield and Reitsma 1991). It might be easy to imagine that multiple quantitative loci that express depending on environment could explain continuous phenotypic plasticity, and there are extensive theoretical studies on continuous plasticity, basically using quantitative genetics models (see below). In contrast, discrete phenotypic plasticity seems rare, indicating that discrete phenotypic plasticity is more difficult to evolve than continuous phenotypic plasticity. Examples include temperature-dependent sex determination (Merchant-Larios and Diaz-Hernandez 2013), male horn of beetle (Nijhout 2003), and a mouth form of nematode (Serobyan et al. 2014). Discrete phenotypic plasticity should arise through development and require a switch that regulates the expression of many genes concurrently, which should be one of the reasons why discrete phenotypic plasticity hardly evolves. In this study, we present a theoretical model aimed at investigating the evolutionary dynamics of such a system and delineating the conditions under which it could evolve. Our approach employs a population genetic framework that integrates the developmental process, where each individual undergoes growth to develop its phenotype, which subsequently becomes subject to selection pressures. Specifically, we concentrate on elucidating the evolutionary trajectory of the gene regulatory network underlying the developmental process, particularly towards the emergence of discrete phenotypic plasticity.

Fig. 1. Illustration of discrete and continuous phenotype plasticity.

While this article specifically focuses on discrete plasticity, the majority of previous theoretical studies on phenotypic plasticity have utilized quantitative genetics models, primarily emphasizing continuous plasticity. These models typically assume the presence of numerous loci with small effects that express themselves differentially depending on environmental conditions. Consequently, the collective effects of these quantitative loci give rise to a continuous distribution of phenotypes, making discrete plasticity highly improbable. Within the framework of quantitative genetics models, adaptive plasticity typically evolves in response to either temporally or spatially fluctuating environments (Via and Lande 1985; Gomulkiewicz and Kirkpatrick 1992; Gavrilets and Scheiner 1993). Numerous factors constraining the evolution of plasticity have been extensively explored, including strong genetic constraints (Via and Lande 1985; Gomulkiewicz and Kirkpatrick 1992), costs associated with plasticity (Van Tienderen 1991, 1997), and the overall unpredictability of environmental conditions (Moran 1992; Gavrilets and Scheiner 1993; De Jong 1999; Tufto 2000; Scheiner and Holt 2012). Additionally, some theories have highlighted the potential role of plasticity in rapidly changing environments (Lande 2009; Chevin and Lande 2010; Reed et al. 2010; Scheiner et al. 2017).

To emphasize discrete plasticity, this study employs a fundamentally different model that integrates the developmental process within each individual. This developmental process is governed by a gene regulatory network, which is believed to play a crucial role in regulating plastic traits in empirical systems (Nijhout 2003; Projecto-Garcia et al. 2017; Ledón-Rettig and Ragsdale 2021). In contrast to previous quantitative genetics models, which rely on statistical correlations between phenotype and environments without considering underlying genetic mechanisms (Schlichting and Pigliucci 1993), our model is more mechanistic-based and plausible. By incorporating a gene regulatory network, we aim to capture the mechanistic basis of phenotypic plasticity, reflecting the biological reality of how genes interact to influence phenotypic expression.

We drew inspiration from the work of Watson et al. (2014), who demonstrated the “memorizing” capability of gene regulatory networks. In their study, the authors simulated the evolution of a gene regulatory network in alternating environments (E1→E2→E1→E2⋯), each with its specific optimal adult phenotype (i.e. expression pattern). As the network evolved, they gradually adapted to both environments, ultimately acquiring the ability to produce either of the two optimum adult phenotypes based on the genotype of the embryo. This resulted in a discrete distribution of phenotypes, where even slight changes in the embryo genotype (by mutation) could lead to vastly different adult phenotypes. In their model, phenotypic variation is clearly discrete. Here, the network’s input, represented by the embryo expression pattern, functions as a switch that determines the resulting adult phenotype. The reason why this ability can be acquired is because the network has experienced the two environments many times in the past, therefore it evolved such that it can adapt to the two environments with minimum changes in its input.

We consider that such a memory ability of gene regulatory networks should be the mechanism behind discrete phenotypic plasticity. If so, a simple modification of the model of Watson et al. (2014) allows to make a realization of phenotypic plasticity with discrete phenotypes. It should be noted that the model of Watson et al. (2014) requires mutation to switch the adult genotype. Here, in order to model phenotypic plasticity as defined, our model is designed such that no mutation is needed to switch the adult genotype. To this end, our gene regulatory network model plugs in a sensor gene, whose expression pattern changes depending on the environment. If we consider the sensor gene mimics hormone receptors, this setting should be reasonable according to empirical evidence: It has been reported that hormone receptors play important roles in gene regulatory systems that express only in a sensitive period in development and work as transcription factors when the corresponding hormone exists in that period (Nijhout 2003; Projecto-Garcia et al. 2017; Ledón-Rettig and Ragsdale 2021). We demonstrate the successful evolution of a sensor gene into a switch, resulting in the emergence of discrete adaptive phenotypes without alterations in genotype. While our approach may appear straightforward, to our knowledge, there are currently no models that precisely elucidate the mechanisms underlying the evolution of discrete phenotypic plasticity through the modulation of gene regulatory networks (as discussed below).

The evolution of gene regulatory networks has been extensively studied in a wide range of literature. A model of gene regulatory network was first developed by A. Wagner (Wagner 1994, 1996) (see also Kauffman 1969). Since then, it has been widely used to investigate how their structural and functional properties change through evolution in static (Wagner 1996; Siegal and Bergman 2002; Azevedo et al. 2006) or fluctuating environments (Crombach and Hogeweg 2008; Draghi and Wagner 2009; Kaneko 2012; Watson et al. 2014). Only a few studies have considered phenotypic plasticity by incorporating sensors into the network (Draghi and Whitlock 2012; Iwasaki et al. 2013; Brun-Usan et al. 2021; Burban et al. 2022). These studies showed that the evolution of phenotypic plasticity generates a correlation of phenotypic variation between traits (Draghi and Whitlock 2012; Brun-Usan et al. 2021) and affects the amount of phenotypic variation observed among populations that encounter novel environments (Iwasaki et al. 2013). Plasticity also affects how gene regulatory networks evolve during domestication processes (Burban et al. 2022). However, these studies considered only continuous phenotypic plasticity (but see Iwasaki et al. 2013). In addition, previous studies have predominantly focused on the outcomes following the establishment of phenotypic plasticity, rather than on the process of its acquisition. Our work stands out by demonstrating how a network can gradually evolve to acquire discrete plasticity, with a crucial role played by the memory ability of the network. We then explore the conditions necessary for a network to achieve such discrete plasticity.

Model

We expand upon the model presented by Watson et al. (2014) by introducing a sensor gene that functions as an environmental detector. Specifically, we examine two distinct environments, denoted as E1 and E2, characterized by a single parameter such as temperature. The sensor gene is active only during a sensitive developmental period and its expression level varies depending on the environmental context. For simplicity, we assume that the sensor gene is active in E1 and inactive in E2. By integrating this sensor gene into the model, our objective is to investigate whether the system can develop phenotypic plasticity, enabling it to produce distinct phenotypes in response to environmental changes, all without altering the genotype. Furthermore, while (Watson et al. 2014) employed a deterministic approach, we enhance realism by incorporating stochasticity, specifically random genetic drift, within the context of finite population genetics.

Our model considers a population of haploid species, in which each individual has three types of genes, G, B, and C. There are n type-G genes and the phenotype of an individual is defined by their expression pattern, which changes during development. gi represents the genotype of the ith type-G gene (i=1,2,3,…,n), which directly determines the expression level of this gene in the embryo (xi(0)=gi). xi changes through development, and let xi(τ) be the expression level of the ith type-G gene at the developmental step τ. We assume that an embryo becomes an adult after τv steps and lives until τ=τl. Natural selection works on the expression pattern of the n type-G genes, x→(τ)=(x1(τ),x2(τ),…,xn(τ)) during adult phase (i.e. τv≤τ≤τl). The types-B and -C genes are involved in this developmental process. bij determines the regulatory effect of jth type-G genes on the expression level of the ith type-G genes. Biologically, the type-B genes can be considered as genomic region(s) that consists of regulatory genes/elements of the type-G genes. In addition, the type-C gene determines how the sensor gene affects the expression of the n type-G genes (ci denotes the effect on the ith type-G gene). Through the developmental process, the expression level of type-G genes changes according to the following formula:

(1) xi(τ+1)=(1−α)xi(τ)+f(∑j=1nbijxj(τ)+ciδ(τ)),

where α is the decay rate of gene expression per developmental step. f is an activation function: we assume that f is a sigmoid function, f(x)=(1+tanh(x))/2, which is a commonly used function to transform any value in (−∞,∞) into the (0,1) interval. In this settings, the minimum and maximum values of xi are 0 and 1/α, respectively. δ(τ) is the expression level of the sensor gene at time step τ. In order to make the effect of the sensor gene minimal so that the developmental process mainly depends on the type-B genes, we assume that the sensor gene expresses only in the first τc steps of the whole developmental process with τl steps (i.e. For τ<τc, δ(τ)=1/α in E1 and δ(τ)=0 in E2). It is important to notice that the absolute expression levels matter only for the type-G genes, on which selection works. The types-B and -C genes determine how the expression level of type-G genes changes over development.

Using this model, we allow all three types of genes to evolve by mutation, selection, and genetic drift, where the environment changes between E1 and E2 at an interval of IT generations. We then ask if the species can display phenotypic plasticity, which is defined as the ability to develop an adult phenotype which is fit to environment E1 if grown in E1 and fit to E2 when grown in E2 from the same genotype (as illustrated in Supplementary Fig. S1). Let X→i=(Xi1,…,Xin) be the optimal expression levels of the n type-G genes in environment Ei (Xij∈[0,1/α]). The fitness of an adult with x→(τ) in environment Ei at developmental step τ is assumed to be given by exp(−|x→(τ)−Xi→|2/(2σ2)). Then, the lifetime fitness is defined as its geometric mean during the adult phase:

(2) wi=[∏τ=τvτlexp(−|x→(τ)−X→i|2/(2σ2))]1/(τl−τv+1)=∏τ=τvτlexp(−∑j=1n[xj(τ)−Xij]22σ2(τl−τv+1)).

The evolution of the three types of genes proceeds by mutation and selection in a Wright–Fisher haploid population of size N. In all individuals, as the initial genotype state, we assume gi=1/(2α),bij=0, and ci=0. Mutation occurs in all three types of genes. μg is the mutation rate per type-G gene per individual per generation, which is constant for all n type-G genes. A mutation changes gi by Δg, where Δg is a random value drawn from a uniform distribution U(−γg,γg). We assume 0≤gi≤1/α, so that gi=1/α if it exceeds 1/α whereas gi=0 if gi is smaller than 0. The mutation rate μb per generation applies to all bij of the type-B genes, and each mutation changes bij by Δb where Δb is drawn from U(−γb,γb). The mutation rate of each ci of the type-C genes is μc per generation, and a mutation changes it by Δc, where Δc is drawn from U(−γc,γc). At each generation, we compute the fitness of all individuals, and the individuals of the next generation are chosen (free recombination is assumed), to which new mutations are introduced.

Through a simulation run, we compute both w1 and w2 every generation, where wi is the fitness of the adult individual if it has grown in environment Ei. The population averages are denoted by w¯1 and w¯2, respectively. If the environment is E1, only w1 is involved in selection and w2 is evolutionary meaningless, and vice versa. Nevertheless, it is interesting to monitor both w1 and w2 through a simulation run to measure the potential ability of the system to fit the two environments.

Results

Overview

To investigate the acquisition of phenotypic plasticity through evolution, we conducted comprehensive simulations utilizing the model outlined earlier. Our primary focus is on scenarios with n=8 (referred to as the full model), although we also explore a simpler model featuring two type-G genes (referred to as the simplified model, as detailed below). As the initial state at T=0 for each simulation run in the full model, we set g=(2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5), c=(0,0,0,0,0,0,0,0), and bij=0, therefore the fitness (w1 and w2) is relatively low. When both w1 and w2 increase and become stably close to 1, we consider that phenotypic plasticity is attained. In each run, we simulated at most 107 generations and terminated when phenotypic plasticity successfully evolved. We considered that phenotypic plasticity had evolved at generation T=Tp if the population mean fitness in both environments (w¯1 and w¯2) was >0.95 in the following 105 generations. If phenotypic plasticity was not acquired until T=107, we considered that plasticity has not evolved successfully. In the following, we consider the parameter values listed in Table 1 unless otherwise specified.

Table 1. Default parameter values in the simulation under the full model with eight type-G genes.

Parameters	Explanation	Default	
N	Population size	1,000	
α	The decay rate of gene expression	0.2	
τc	The developmental step when the sensor gene expression ends	10	
τv	The developmental step when the expressions start to affect fitness	21	
τl	The developmental step when the expressions end to affect fitness	40	
σ	Inverse of the strength of selection	5	
μi	Mutation rate per gene for the type-i gene	10−5	
γi	Maximum effect size of a new mutation for the type-i gene	0.1	
IT	The number of generations between environmental switch	—	
X→1	Optimal expression pattern in E1	(0, 0, 5, 5, 5, 5, 0, 0)	
X→2	Optimal expression pattern in E2	(0, 0, 5, 5, 0, 0, 5, 5)	

Figure 2 summarizes the simulation results for various pairs of mutation rates and IT under the full model with eight type-G genes. We performed 100 independent simulation runs for each parameter set, from which the proportion of runs that successfully acquired plasticity, Sp, and Tp in the successful runs were recorded. Figure 2 demonstrates that phenotypic plasticity successfully evolves for almost all runs. The waiting time for the acquisition of plasticity (Tp) tends to be short at intermediate IT (IT=100,1,000). It can be explained by the difference of the evolutionary dynamic between small and large IT. In the cases of a very small IT, the network should perform a challenging task to adapt to both environments almost simultaneously. For a large IT, the network can adapt to the environments in turn, perhaps making the adaptation process straightforward. However, there is a drawdown for a large IT. When IT is very long, signifying an extended period during which the system adapts to one environment, it becomes challenging to maintain the ability to produce a phenotype suited for the other environment. However, despite this challenge, plasticity eventually evolves after a substantial number of environmental changes, albeit over a prolonged duration. Our findings illustrate that given a sufficiently long time, even with a large IT, such as IT≥104 generations, the system possesses significant potential to retain past adaptations. This memory capability enables the system to appropriately switch phenotypes based on the prevailing environment.

Fig. 2. Summary of simulation results in the full model. The effects of the three mutation rates a) μg, b) μb, and c) μc) on the proportion of successful runs (Sp) and waiting time (Tp) are shown for various IT. In each panel, Sp is shown in the upper bar plot while the distribution of Tp in successful runs is shown in the lower violin plot. For each μi, five μi values are examined. All other parameters are set to their default values (see Table 1). As a control, the middle column (μi=10−5) presents the results for the default parameters, so that the results are identical in (a), (b), and (c).

Although we have provided a brief overview of the simulation results, understanding the impact of each parameter requires examining how the system evolves and ultimately acquires plasticity. To facilitate this understanding, we will first examine the evolutionary process in typical simulation runs. Given the complexity of the dynamics involved in the full model with eight type-G genes, we will begin with a simplified model featuring two type-G genes for elucidation purposes. Subsequently, we will verify that similar dynamics are also observed in the full model. Following a description of the evolutionary mechanism, we will revisit Fig. 2 to discuss the influence of each parameter on the evolution of plasticity.

Evolutionary dynamics in the simplified model

To understand how the system acquires plasticity, we performed simulations under a very simplified setting with only two type-G genes (n=2). The optimum expression levels of the two type-G genes are X→1=(5,0) in E1 and X→2=(0,5) in E2. We assumed that a higher mutation rate of μi=10−4 and stronger selection σ=5/2 to account the reduced number of genes compared with the full model. All other parameters are same as those in Table 1.

Figures 3 and 4 detail typical behaviors of the evolutionary process towards plasticity. The process is quantitatively quite different when IT is very small and when IT is very large, therefore we consider two cases, IT=1 and 10,000. The adaptation process to multiple environments in each run can be easily understood by examining the trajectories of the two fitness values, w¯1 and w¯2. We first consider the case of IT=1, where the environment changes every generation and the system has to adapt to both environments almost simultaneously. Figure 3 shows the evolutionary trajectory in a represented run, in which, after w¯1 and w¯2 stay for a while at w¯1∼w¯2∼0.6, they jump up and reach a plateau at w¯1∼w¯2∼1 (see Fig. 3a).

Fig. 3. Evolution of phenotypic plasticity in the simplified model with two type-G genes when IT=1. The result of a single representative run is illustrated. a) Temporal change of fitness through evolution. The red and blue lines are the fitness in environment E1 and that in E2, respectively (i.e. w¯1 and w¯2). b) Phase portraits of the development of a single individual at five time points (T=0,80,000,90,000,100,000, and 400,000) shown by the dashed vertical lines in (a). (x1(τ), x2(τ)) starts at × at embryo and moves along the red and blue arrows in E1 and E2, respectively. The star represents (x1(τc), x2(τc)) when the activation of the switch ends in E1. The default parameter values were used (see Table 1) except for μi=10−4, σ=5/2, X→1=(5,0), and X→2=(0,5).

Fig. 4. Evolution of phenotypic plasticity in the simplified model with two type-G genes when IT=10,000. The result of a single representative run is illustrated. a) Temporal change of fitness through evolution. The red and blue lines are the fitness in environment E1 and that in E2, respectively (i.e. w¯1 and w¯2). b) Phase portraits of the development of a single individual at five time points (T=0,50,000,60,000,122,000, and 400,000) shown by the dashed vertical lines in (a). (x1(τ), x2(τ)) starts at × at embryo, and moves along the red and blue arrows in E1 and E2, respectively. The star represents (x1(τc), x2(τc)) when the activation of the switch ends in E1. The default parameter values were used (see Table 1) except for μi=10−4, σ=5/2, X→1=(5,0), and X→2=(0,5).

Figure 3b shows how the expression levels of the two type-G genes (x1 and x2) changes along development in each individual at five time points. At T=0, we set g=(2.5,2.5), c=(0,0), and bij=0, therefore at embryo, (x1(0), x2(0)) locates at the center (2.5, 2.5). As this is a stable equilibrium in the developmental dynamics defined by the type-B genes (bij), (x1(τ), x2(τ)) does not change along development, resulting in (x1(τ), x2(τ))=(2.5, 2.5) for any τ.

At T=80,000, with w¯1∼w¯2∼0.6 (time point ▽T1a in Fig. 3), both w¯1 and w¯2 begin to increase as bij, ci, and gi evolve, bringing the adult phenotype closer to its optimum in both environments. At this stage, there is only one stable equilibrium close to the center in the phase portrait (depicted by the purple circle), indicating that individuals develop towards this equilibrium regardless of the environment. However, the difference between the two environments lies in the developmental trajectory: in E2, (x1(τ), x2(τ)) moves directly towards the purple circle (as indicated by the blue arrow), while in E1, the evolved switch causes (x1(τ), x2(τ)) to initially jump to the star before eventually reaching the equilibrium. A similar pattern is observed at T=90,000 (time point ▽ T1b), except that the limited number of developmental steps prevents (x1(τl), x2(τl)) from reaching the equilibrium. During this phase, the network exhibits only one equilibrium, indicating a phenotype that is moderately adapted to both environments. It appears that the network is still struggling with how to fully adapt to the alternating environments.

A drastic change arises at T=100,000 (time point ▽ T1c), where both w¯1 and w¯2 are close to 1. This means that plasticity is almost achieved, and the reason why a single genotype can design two distinct phenotypes is as follows. The key is that the network has evolved such that Equation (1) defined by the type-B genes has two stable equilibria, one close to the top-left corner (blue circle) and the other is the bottom-right corner (red circle), which are obviously the two optimum phenotypes in the two environments. A dashed line is introduced such that (x1(τ), x2(τ)) transitions towards the blue circle when above the line and towards the red circle when below it (i.e. basin boundary). Consequently, if an individual develops in E2, (x1(τ), x2(τ)) directly converges to the equilibrium near the optimal blue circle. Conversely, in E1, where the switch is activated, (x1(τ), x2(τ)) initially crosses the line towards the star, subsequently transitioning automatically to the red circle. In this scenario, plasticity is deemed to be nearly acquired. As the system evolves further, the two equilibria converge near (0,5) and (5,0), where both w¯1 and w¯2 approach 1 (T=400,000). While the expression levels of the embryo (x1(0), x2(0)) have shifted slightly from the center, their contribution to the evolution of plasticity diminishes.

Figure 4 demonstrates how phenotypic plasticity evolves with a very large IT. Starting at the same initial state as in Fig. 3, the system immediately reaches a phase where w¯1 and w¯2 oscillate dramatically (Fig. 4a). For instance, in E1 at time point ▽ T1a, the equilibrium approaches (5,0), indicating adaptation solely to E1. Consequently, the developmental outcome is uniform across environments, with development consistently progressing to the single equilibrium (purple circle). After an environment change, the system tries to adapt to E2 and the equilibrium moves to (0,5) at time point ▽ T1b.

Repeated adaptation to the two environments leads to a qualitative change in the phase portrait. For instance, consider time point ▽ T2, at which the system is in E1. Although the developmental outcome remains consistent across environments, we observe an additional equilibrium near the optimal expression for E2 (i.e. the top-left corner). This new equilibrium is a consequence of the system’s prior adaptation to E2, which persists even in the current phase within E1, albeit unused. This scenario echoes the concept of “memorization” of past adaptations as described by Watson et al. (2014). By time point ▽ T3, the network has acquired phenotypic plasticity through the successful evolution of the switch function of the type-C gene. This allows (x1(τ), x2(τ)) to traverse the dashed line when the switch is activated. Consequently, the two equilibria are well-adapted, situated very close to (0,5) and (5,0).

Thus, our simulation under the simplified model effectively demonstrates the fundamental mechanism underlying the acquisition of discrete phenotypic plasticity. Two key conditions appear necessary for achieving phenotypic plasticity: (1) The gene regulatory network must possess a memory of past adaptations, manifesting as two stable equilibria corresponding to the two optimal expression patterns. (2) The type-C gene must facilitate the traversal of the developmental phenotype (expression) across the boundary of the basin during development. With these conditions in place, our model inherently generates a system in which only two distinct phenotypes emerge. However, there may be exceptions, such as instances where some individuals fail to reach an equilibrium in the developmental process due to slow convergence. In these cases, intermediate phenotypes may be observed, potentially leading to a nondiscrete distribution of phenotypes.

It would be intriguing to explore the scenario where an evolved individual with acquired plasticity is exposed to an intermediate environment between E1 and E2. In such a situation, we would expect the switch to be partially activated. Nevertheless, the system would still produce two phenotypes as long as the developmental process occurs rapidly enough for all individuals to develop optimal phenotypes. Intermediate phenotypes could arise if the partially activated switch impedes (x1(τ), x2(τ)) from effectively crossing the basin boundary (i.e. if the point represented by the star is located near the boundary), preventing it from reaching the optimum before reaching adulthood. To explore these possibilities further, we will return to the full model with eight type-G genes to consider more realistic scenarios.

Cases of fast environment transition in the full model

Following Figs. 3 and 4, we trace the adaptation process to multiple environments in each run by examining the trajectories of the two fitness values, w¯1 and w¯2. We first consider the case of IT=1, where the environment changes every generation and the system has to adapt to both environments almost simultaneously. The overall pattern is very similar to that shown in Fig. 3. Figure 5 shows the evolutionary trajectory in a represented run, where both w¯1 and w¯2 increase simultaneously in the earlier phases, and after staying for a while at w¯1∼w¯2∼0.6 they increase and achieve w¯1∼w¯2∼1 (see Fig. 5a). Since it is no longer feasible to make a phase portrait in the full model, we below see the evolutionary dynamics focusing on the genotype and adult phenotype.

Fig. 5. Evolution of phenotypic plasticity in the full model when IT=1. The result of a single representative run is illustrated. a) Temporal change of fitness through evolution. The red and blue lines are the fitness in environment E1 and that in E2, respectively (i.e. w¯1 and w¯2). b) Genotype and phenotype of an adult in each environment at three time points (T=0,2×105, and 2×106) shown by the dashed vertical lines in (a). Phenotypes at the middle of the adult phase τ=30 are shown. The default parameter values were used (see Table 1).

In this run, at T=0, we set g=(2.5,2.5,2.5,2.5,2.5,2.5,2.5,2.5), c=(0,0,0,0,0,0,0,0) and bij=0. With this setting, the adult phenotype at T=0 is x→(τ)=g whichever environment it is grown in (all genes are in brown in Fig. 5b), which are quite different from the two optimum phenotypes so that the fitness is quite low (w¯1∼w¯2∼0.3) (time point ▽ T0 in Fig. 5). w¯1 and w¯2 then start increasing along the evolution of bij, ci and gi to make adult phenotype close to its optimum in both environments. At T=2×105 with w¯1∼w¯2∼0.6 (time point ▽T1 in Fig. 5), the expression patterns of the first four type-G genes (genes 1, 2, 3, and 4) are already very close to the optimums, which is easy to understand because their optimum expression patterns are shared by the two environments. In other words, the evolution of these genes are rather straightforward with a single goal (i.e. optimum expression), which does not apply to the other four genes with two distinct goals. At time point ▽T1, genes 5, 6, 7, and 8 have not adapted to either of the two environments, exhibiting similar adult expression levels in the two environments. In other words, the network can produce only one kind of phenotype that are moderately adapted to the two environments, consistent with time point ▽ T1a in Fig. 3.

The network eventually acquires the ability to produce two optimum phenotypes almost completely around time T=4×105 (see Fig. 5a). At time T=2×106 (time point ▽T2 in Fig. 5), genes 5, 6, 7, and 8 have distinct expression levels depending on the environment, resulting in two distinct phenotypes. The sensor gene (type-C gene) has evolved as an initial trigger for developmental plasticity by enhancing type-G genes 5 and 6 and suppressing type-G genes 7 and 8 in E1. In E2, high expression of type-G genes 7 and 8 and very low expression of type-G genes 5 and 6 are produced without the expression of the sensor gene (directly through the regulation network of type-B genes itself). The relative contribution of type-G genes to the acquisition of plasticity seems low.

Similar patterns were observed in almost all runs (see Supplementary Fig. S2 for other runs). It can be summarized that the adaptation process generally involves two steps when IT is small. In the first step, the network evolves such that it produces a single kind of phenotype (whichever the sensor gene is active or inactive) that is equally fit to the two environments, but the fitness is fairly low. Then, the network struggles for a while to evolve to be able to produce two optimum phenotypes. During this process, we do not observe a marked increase in fitness because fitting to one environment likely has a negative effect on fitting to the other environment. In the meantime, the network still keeps evolving, and at some point, it becomes able to produce two optimum phenotypes when both w¯1 and w¯2 dramatically increase to almost 1 (second step). We can consider that, at this stage, the network has “memorized” the two environments, obtaining the ability to produce two adaptive phenotypes with a help of the switch function by the sensor gene.

Cases of slow environment transition in the full model

Next, we focus on how phenotypic plasticity evolves with a very large IT. Figure 6 shows the result of a representative run when IT is very large (IT=10,000). The pattern is quite different from that for a small IT in the previous section. Starting at the same initial state as in Fig. 5, the system immediately reaches a phase where w¯1 and w¯2 dramatically fluctuate (Fig. 6). In the earlier stages in this phase, the fitness increases gradually after each environment change, and then the fitness’ recovery becomes faster over time. As is clearly seen in the close-up of Fig. 6a, after an environment change, for instance from E1 to E2 (e.g. at time point ▽ T1a), w¯2 suddenly increases to ∼1 whereas w¯1 decreases from ∼1 (i.e. losing the ability to adapt to E1), and vice versa (e.g. at time point ▽ T1b). The observed immediate adaptation suggests that the system has evolved the ability to adapt to the environment change very quickly through a few mutations after experiencing adaptation to the same environment, consistent with the result of Crombach and Hogeweg (2008). This means that, during this phase with fluctuating w¯1 and w¯2, the network gradually “memorizes” how to adapt to each of the two environments. At the moment, the network can produce two phenotypes by a minimum number of mutations, but the sensor gene does not work well as a switch (i.e. phenotypic plasticity has not evolved yet).

Fig. 6. Evolution of phenotypic plasticity in the full model when IT=10,000. The result of a single representative run is illustrated. a) Temporal change of fitness through evolution. The red and blue lines are the fitness in environment E1 and that in E2, respectively (i.e. w¯1 and w¯2). b) Genotype and phenotype of an adult in each environment at four time points (T=0,2.9×105,3×105, and 2×106) shown by the three dashed vertical lines in (a). Phenotypes at the middle of the adult phase τ=30 are shown. The default parameter values were used (see Table 1).

The middle two panels in Fig. 6b show the expression patterns at time points ▽ T1a and ▽ T1b, where the average of w¯1 and w¯2 is roughly 0.6. At time point ▽ T1a (the end of E1), the system is fully adapted to E1, but not yet to E2. The expression patterns of the type-G genes in E1 and E2 are very similar to each other and fit to only E1. After IT=10,000 generations at time point ▽ T1b (the end of E2), the expression patterns of the type-G genes in E1 and E2 are again very similar to each other and fit to only E2. This phase is comparable with time points ▽ T1a and ▽ T1b in the simplified model (Fig. 4). It is interesting to note that the clear difference in the expression patterns of type-G genes between ▽T1a and ▽T1b is caused by a very little change in bij, ci, and gi, and no notable difference is observed. The network still needs to evolve such that it can switch the two phenotypes even more easily, ideally by the sensor gene alone (without changing the network itself by mutation). The phase of alternate fluctuation of w¯1 and w¯2 continues for quite long, then the system suddenly acquires the ability to fit to both environments (w¯1∼w¯2∼1) around T=8×105. Even after this event, we can observe a drastic drop of fitness, but it recovers very quickly.

The bottom panel of Fig. 6b shows the expression patterns at time point ▽ T2, where the system successfully has adapted to the two distinct optimums X→i without changing bij, ci and gi. The types-B and -C genes have contributed to this acquisition of plasticity, which is seen in well activated bij and ci for i,j≥5. The relative contribution of type-G genes to the acquisition of plasticity seems low.

This typical pattern is observed in almost all runs using this parameter set (see Supplementary Fig. S3). The pattern is quite different from the case of a small IT where the system attempts to adapt to the two environments simultaneously, indicating IT affects how a network acquires plasticity. It seems that, while repeatedly adapting to one of the two environments, the network gradually “memorizes” how to adapt, and becomes possible to adapt to environmental changes with minimum modifications by mutation. Then, the sensor gene eventually evolves as an effective switch, which offers an ability to adapt the two environment with no mutational changes in the system.

Thus far, we showed two typical patterns for small and large IT. We found that, in cases with intermediate IT, the pattern is a mixture of the two extreme cases. See Supplementary SI for details.

Effect of mutation parameters

Based on the detailed understanding of the evolutionary dynamics toward discrete plasticity (see above), we explore the effect of mutation parameters in the full model with eight genes. In Fig. 2, we summarize how each of the three mutation rates (μg, μb, μc) affects Sp and Tp. In general, they have small effects on Sp, while Tp is quite much affected by the mutation rates. μc is the mutation rate at the switch gene. Since quick changes of the switch gene enhances the ability to react environmental changes and promote the acquisition of plasticity, Tp decreases as μc increases, irrespective of IT. The mutation rate of the type-B genes (μb) has a more complicated effect: μb decreases Tp when IT is small, while it slightly increases Tp when IT is large. This pattern is explained by the difference in the evolutionary dynamics between the small and the large IT cases (see Figs. 5 and 6). When IT is small, the network needs to adapt to the two environments almost simultaneously, which requires joint evolution of the type-B and the type-C genes. Since the faster evolution of the type-B genes makes the acquisition of this ability easier (as well as type-C-gene), μb has a negative correlation with Tp. In contrast, when IT is large, the gene regulatory network adapts to each environment one by one. This process can be rephrased as follows. When the network tries to adapt to one environment, it has to modify the network that was optimized to the other environment. That is, learning (adapting to) one environment is coupled with reducing the fitness to the other (i.e. forgetting the past adaptation). Mutations in the type-B genes not only promotes adaptation to the current environment but also losing adaptation to the previous environment. Therefore, when IT is large, μb has a slightly positive correlation with Tp. The effect of μg seems small, suggesting that gi is relatively less important in the evolution of plasticity.

We also examined the effect of mutational impacts (γg, γb, γc) while assuming μg=μb=μc=10−5 (Supplementary Fig. S4). The basic pattern of the effect of γi is very similar to that of μi, although changes in γi have a greater impact on Tp and also have some impact on Sp. As γc increases, plasticity is more likely to evolve in a short time. A large γb tends to decrease Tp when IT is small, while increases Tp when IT is large. Sp becomes small in the parameter space where Tp is very long (i.e. small γc, and large γb in large IT). No clear effect of γg is observed, consistent with the absence of the effect of μg.

Discreteness of phenotype

Our model thus far considered two distinct environments, E1 and E2, each of which has its optimum phenotype, and a successful adaptation to both environments is required to evolve plasticity with discrete phenotypes. It was assumed that the sensor gene is active in E1 and inactive in E2, or the expression level is 1/α and 0 in E1 and E2, respectively. To check whether the evolved plasticity is discrete or continuous, it would be intriguing what kind of phenotype arises if the environment is somewhere between E1 and E2. We here consider intermediate environments between E1 and E2, which can be realized by assuming intermediate levels of expression of the sensor gene between 0 and 1/α.

With this simple assumption, we re-examined the simulation results. Each simulation run produced an evolved gene regulatory network (types-G, B, C genes), which can be used to reproduce adult phenotypes under any environment. In practice, after each run, for each individual, a large number of its embryos were grown in intermediate environments. We then judged that the distribution is discrete if only two optimum phenotypes arise in almost all range of the environment. Figure 7a shows a typical example of an individual with discrete plasticity. Focus on the expression patterns of genes 5, 6, 7, and 8 that distinguish the two optimum phenotypes, X→1 and X→2. When the environment is changed from 0 to 1 (i.e. from E1 to E2), we observe that the expression levels of these genes abruptly change between 0 and 5 around when the environmental parameter is 0.75, indicating that we can consider that only two types of phenotypes arise in almost all range of the interval. This pattern is also well-characterized by the distribution of w1 and w2: When the distribution is discrete, in almost all range, either w1 or w2 is nearly 1, with only a small range where both w1 are w2 are significantly smaller than 1. In Fig. 7c, a typical pattern of the cases of continuous phenotypes is shown. When the environment is transitioned from 0 to 1, the expression levels of genes 5, 6, 7, and 8 change gradually. As a consequence, w1 and w2 change gradually and there is a fairly large range of the environment that produces “intermediate phenotypes” (i.e. phenotypes that does not fit to either E1 or E2).

Fig. 7. Discreteness of phenotype distribution in intermediate environment. a, b) An example case with a highly discrete distribution. Expression patterns of the 8 type-G genes a) and fitness values in E1 and E2 b) are plotted. Default parameters and IT=1 were assumed. c, d) An example case with a less discrete distribution. Default parameters except for μc=10−6 and IT=1 were assumed.

Based on this logic, we examined if the distribution of phenotype is discrete by using all successful runs in Fig. 2. We defined that intermediate phenotypes are those with both w1 and w2 less than 0.95. We checked the phenotypes of all individuals in the entire range of the environment and considered that the plasticity is discrete if intermediate phenotypes are less than 5%. Figure 8 shows the proportion of runs with a discrete phenotype distribution in each parameter set, where one of the three mutation rate parameters (μg,μb,μc) was changed, while the other two were fixed. It was found that the proportion of discrete plasticity varies depending on the parameter sets. The proportion is large when IT is large, suggesting that the discrete plasticity is more likely to evolve when the network learn the two phenotypes one by one, rather than in a simultaneous manner. The proportion tends to be large when μb and μc are large. This pattern reflects that the well-evolved type-B and type-C genes are needed to make discrete plasticity so that the convergence to each equilibrium is sufficiently fast. The effect of μg is not well observed, consistent with the results in the previous section.

Fig. 8. Proportion of simulation runs with a discrete phenotypic distribution. μg a), μb b), and μc c) were changed. Other parameters are the default values.

Discussion

This study illustrates how discrete phenotypic plasticity could emerge through the evolutionary dynamics of a gene regulatory network. Our simulations in the full model with eight type-G genes showed that phenotypic plasticity could successfully evolve almost all parameter sets investigated (Fig. 2, Supplementary Fig. S4), indicating a marked ability of gene regulatory network to produce multiple phenotypes properly, but it takes a very long time. The simplified model with two type-G genes helped to understand the mechanism behind the acquisition of discrete plasticity (Figs. 3 and 4). (1) The most important aspect is the evolution of the phase portrait of the developmental process, which is orchestrated by the type-B genes. These genes must establish two stable equilibria, each corresponding to the optimal expression pattern under one of the two environments. This ensures that development can progress towards one of the two optima, effectively “memorizing” the adaptation to the environments within the gene regulatory network. (2) The type-C gene functions as a “switch” to enable the phenotype (expression) to cross the basin boundary during the developmental process. This dictates which optimum the development will proceed towards, effectively directing the phenotype without altering the type-B genes. These two factors represent the minimum conditions necessary for the evolution of discrete phenotypic plasticity through the evolution of a gene regulatory network. Additionally, for the phenotype distribution to become more discrete, the developmental process must occur rapidly enough for all individuals to reach their optimum phenotype. Failure to meet this condition may result in some individuals exhibiting intermediate phenotypes. This scenario is particularly relevant when intermediate environments are present, leading to a nondiscrete phenotype distribution along the environmental gradient (see Figs. 7 and 8). This reasoning is consistent with our simulations conducted under the full model with eight type-G genes.

There appear to be two primary pathways in the evolutionary process leading to the acquisition of plasticity, as illustrated in Figs. 5 and 6 (see also Figs. 3 and 4). When the environment changes frequently (i.e. small IT, as shown in Fig. 5), the network initially evolves to produce a phenotype with maximum “average” fitness (i.e. w1w2). Subsequently, both the ability to generate multiple phenotypes and the proper switching mechanism are acquired at a certain point in time. In such scenarios, the co-evolution of the type-B and type-C genes is crucial, and higher mutation rates at these genes accelerate the evolution of plasticity (Fig. 2). In contrast, when IT is large (Fig. 6), the network undergoes repeated adaptations to each environment and initially gains the ability to switch between the two phenotypes with only a few mutations (i.e. memorizing the two environments within a single network). Subsequently, the trigger for phenotype switching shifts from mutations to the sensor gene. Importantly, the network retains a robust memory capability, allowing it to remember the adaptive phenotype to past environments even after evolving for thousands of generations in different environments. In scenarios with slow environmental transitions, a high mutation rate at the type-C genes facilitates the evolution of plasticity by aiding in the acquisition of the proper switching mechanism. Conversely, a high mutation rate at the type-B genes slightly hinders the evolution of plasticity because a faster evolution of the network removes remnants of past adaptations, although it does promote the memory of the two optimum patterns (Fig. 2). Regardless of the specific conditions, we demonstrate that discrete phenotypic plasticity could eventually evolve across wide ranges of parameter sets, albeit requiring a significant amount of time, particularly in cases where IT is large.

Our study successfully demonstrated the evolution of discrete phenotypic plasticity, a feat that has proven challenging in previous models reliant on quantitative genetics, which tend to produce continuous distributions of phenotype (see the Introduction). Verbal arguments have suggested that considering the developmental process is crucial for understanding discrete plasticity, as the action of regulatory switches is necessary to establish such discrete traits (Schlichting and Pigliucci 1993; Via et al. 1995), yet mechanistic-based models of discrete plasticity have been relatively understudied. While some recent studies have explored the evolution of discrete plasticity (Svennungsen et al. 2011; Chevin and Lande 2013), these models have largely remained phenotype-based and have not provided insights into the underlying genetic mechanisms. Our theoretical framework offers a valuable tool for investigating the genetic properties involved in the evolution of discrete phenotypic plasticity, particularly in terms of gene expression, which is expected to differ both quantitatively and qualitatively from continuous plasticity.

This work demonstrates that discrete plasticity is more likely to emerge when mutation rates of the types-B and C genes are relatively high and environmental fluctuations are less frequent. Under these conditions, the convergence to stable equilibria in the developmental process tends to occur reasonably quickly, allowing the phenotype to reach one of the equilibria before adulthood even in intermediate environments, where the sensor gene may be only partially activated (see Fig. 8). This finding may shed light on the relatively rare empirical observations of discrete phenotypic plasticity, as environmental fluctuations in nature tend to occur over shorter time scales, such as seasonal changes. It is important to note that while our model assumes only two environments alternating with equal intervals, environmental changes in nature are much more complex. Given this complexity, the acquisition of discrete plasticity could be even more challenging in natural systems.

Our model considers a gene regulatory network, in which gene members regulate one another, that is, they behave as if they are transcription factors. In practice, transcription factors play the central role in development, and it is known that those genes likely self-regulate so that their expressions should be strongly enhanced or suppressed (Osborn et al. 2011; Bouldin et al. 2015). Provided this, our model would be reasonable in that we set the expression levels to be either 0 or 5 at the optimum states of X→. If we instead assumed intermediate optimum values (e.g. X→1=(1,1,4,4,4,4,1,1) and X→2=(1,1,4,4,1,1,4,4), we found that Sp was reduced and that the evolution of the discrete plasticity is less likely (see Supplementary Figs. S5, S6), suggesting that, in such situations, the evolution of the type-B genes that encompass fast convergence to equilibria in the developmental process is difficult.

In conclusion, our study illustrates that discrete phenotypic plasticity can evolve through the modulation of a regulatory network influenced by environmental cues sensed by a dedicated sensor gene. This modeling approach aligns well with empirical evidence indicating the involvement of hormones in shaping plastic traits during development. In many instances, environmental stimuli impact hormone secretion, while various hormone receptors, acting as transcription factors, regulate the developmental switch governing plastic traits based on hormone concentration levels during critical developmental periods (Nijhout 2003; Projecto-Garcia et al. 2017; Ledón-Rettig and Ragsdale 2021). This work underscores the significance of environmental cues in shaping gene regulatory networks, ultimately leading to the emergence of discrete phenotypic variations associated with developmental plasticity.

Supplementary Material

jkae144_Supplementary_Data

Acknowledgments

We thank two anonymous reviewers for helpful comments. The computing resource was provided by Human Genome Center (the University of Tokyo).

Data availability

Codes used for simulations are available at https://github.com/TSakamoto-evo/plasticity.

Supplemental material available at G3 online.

Funding

This work was partly supported by grants from the Graduate University for Advanced Studies, SOKENDAI, and Japan Society for the Promotion of Science (JSPS) to H.I. (JSPS KAKENHI Grant Numbers JP23H02521, JP23H04736) and T.S. (JSPS KAKENHI Grant Numbers JP20J21760, JP23KJ2158).
==== Refs
Literature cited

Azevedo RB , LohausR, SrinivasanS, DangKK, BurchCL. 2006. Sexual reproduction selects for robustness and negative epistasis in artificial gene networks. Nature. 440 :87–90. doi:10.1038/nature04488.16511495
Bouldin CM , ManningAJ, PengYH, Farr IIIGH, HungKL, DongA, KimelmanD. 2015. Wnt signaling and tbx16 form a bistable switch to commit bipotential progenitors to mesoderm. Development. 142 :2499–2507. doi:10.1242/dev.124024.26062939
Brakefield PM , ReitsmaN. 1991. Phenotypic plasticity, seasonal climate and the population biology of bicyclus butterflies(satyridae) in Malawi. Ecol Entomol. 16 :291–303. doi:10.1111/een.1991.16.issue-3.
Brun-Usan M , RagoA, ThiesC, UllerT, WatsonRA. 2021. Development and selective grain make plasticity ‘take the lead’ in adaptive evolution. BMC Ecol Evol. 21 :205. doi:10.1186/s12862-021-01936-0.34800979
Burban E , TenaillonMI, Le RouzicA. 2022. Gene network simulations provide testable predictions for the molecular domestication syndrome. Genetics. 220 :iyab214. doi:10.1093/genetics/iyab214.34849852
Chevin L , LandeR. 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.19863583
Chevin LM , LandeR. 2013. Evolution of discrete phenotypes from continuous norms of reaction. Am Nat. 182 :13–27. doi:10.1086/670613.23778223
Crombach A , HogewegP. 2008. Evolution of evolvability in gene regulatory networks. PLoS Comput Biol. 4 :e1000112. doi:10.1371/journal.pcbi.1000112.18617989
De Jong G . 1999. Unpredictable selection in a structured population leads to local genetic differentiation in evolved reaction norms. J Evol Biol. 12 :839–851. doi:10.1046/j.1420-9101.1999.00118.x.
Draghi J , WagnerGP. 2009. The evolutionary dynamics of evolvability in a gene network model. J Evol Biol. 22 :599–611. doi:10.1111/jeb.2009.22.issue-3.19170816
Draghi JA , WhitlockMC. 2012. Phenotypic plasticity facilitates mutational variance, genetic variance, and evolvability along the major axis of environmental variation. Evolution. 66 :2891–2902. doi:10.1111/evo.2012.66.issue-9.22946810
Gavrilets S , ScheinerSM. 1993. The genetics of phenotypic plasticity V. Evolution of reaction norm shape. J Evol Biol. 6 :31–48. doi:10.1046/j.1420-9101.1993.6010031.x.
Gomulkiewicz R , KirkpatrickM. 1992. Quantitative genetics and the evolution of reaction norms. Evolution. 46 :390–411. doi:10.2307/2409860.28564020
Huey RB , KingsolverJG. 1989. Evolution of thermal sensitivity of ectotherm performance. Trends Ecol Evol. 4 :131–135. doi:10.1016/0169-5347(89)90211-5.21227334
Iwasaki WM , TsudaME, KawataM. 2013. Genetic and environmental factors affecting cryptic variations in gene regulatory networks. BMC Evol Biol. 13 :91. doi:10.1186/1471-2148-13-91.23622056
Kaneko K . 2012. Evolution of robustness and plasticity under environmental fluctuation: formulation in terms of phenotypic variances. J Stat Phys. 148 :687–705. doi:10.1007/s10955-012-0563-1.
Kauffman SA . 1969. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 22 :437–467. doi:10.1016/0022-5193(69)90015-0.5803332
Lande R . 2009. Adaptation to an extraordinary environment by evolution of phenotypic plasticity and genetic assimilation. J Evol Biol. 22 :1435–1446. doi:10.1111/jeb.2009.22.issue-7.19467134
Ledón-Rettig CC , RagsdaleEJ. 2021. Physiological mechanisms and the evolution of plasticity. In: Pfenig DW, editor. Phenotypic Plasticity & Evolution. Boca Raton (FL): CRC Press. p. 113–137.
Merchant-Larios H , Diaz-HernandezV. 2013. Environmental sex determination mechanisms in reptiles. Sex Dev. 7 :95–103. doi:10.1159/000341936.22948613
Moran NA . 1992. The evolutionary maintenance of alternative phenotypes. Am Nat. 139 :971–989. doi:10.1086/285369.
Nijhout HF . 2003. Development and evolution of adaptive polyphenisms. Evol Dev. 5 :9–18. doi:10.1046/j.1525-142X.2003.03003.x.12492404
Osborn DP , LiK, HinitsY, HughesSM. 2011. Cdkn1c drives muscle differentiation through a positive feedback loop with myod. Dev Biol. 350 :464–475. doi:10.1016/j.ydbio.2010.12.010.21147088
Projecto-Garcia J , BiddleJF, RagsdaleEJ. 2017. Decoding the architecture and origins of mechanisms for developmental polyphenism. Curr Opin Genet Dev. 47 :1–8. doi:10.1016/j.gde.2017.07.015.28810163
Reed TE , WaplesRS, SchindlerDE, HardJJ, KinnisonMT. 2010. Phenotypic plasticity and population viability: the importance of environmental predictability. Proc R Soc B Biol Sci. 277 :3391–3400. doi:10.1098/rspb.2010.0771.
Scheiner SM , BarfieldM, HoltRD. 2017. The genetics of phenotypic plasticity XV. Genetic assimilation, the Baldwin effect, and evolutionary rescue. Ecol Evol. 7 :8788–8803. doi:10.1002/ece3.2017.7.issue-21.29152178
Scheiner SM , HoltRD. 2012. The genetics of phenotypic plasticity. X. Variation versus uncertainty. Ecol Evol. 2 :751–767. doi:10.1002/ece3.217.22837824
Schlichting CD , PigliucciM. 1993. Control of phenotypic plasticity via regulatory genes. Am Nat. 142 :366–370. doi:10.1086/285543.19425982
Serobyan V , RagsdaleEJ, SommerRJ. 2014. Adaptive value of a predatory mouth-form in a dimorphic nematode. Proc R Soc B Biol Sci. 281 :20141334. doi:10.1098/rspb.2014.1334.
Siegal ML , BergmanA. 2002. Waddington’s canalization revisited: developmental stability and evolution. Proc Natl Acad Sci USA. 99 :10528–10532. doi:10.1073/pnas.102303999.12082173
Svennungsen TO , HolenØH, LeimarO. 2011. Inducible defenses: continuous reaction norms or threshold traits? Am Nat. 178 :397–410.21828995
Tufto J . 2000. The evolution of plasticity and nonplastic spatial and temporal adaptations in the presence of imperfect environmental cues. Am Nat. 156 :121–130. doi:10.1086/303381.10856196
Van Tienderen PH . 1991. Evolution of generalists and specialists in spatially heterogeneous environments. Evolution. 45 :1317–1331. doi:10.1111/j.1558-5646.1991.tb02638.x.28563821
Van Tienderen PH . 1997. Generalists, specialists, and the evolution of phenotypic plasticity in sympatric populations of distinct species. Evolution. 51 :1372–1380. doi:10.1111/evo.1997.51.issue-5.28568610
Via S , GomulkiewiczR, De JongG, ScheinerSM, SchlichtingCD, Van TienderenPH. 1995. Adaptive phenotypic plasticity: consensus and controversy. Trends Ecol Evol. 10 :212–217. doi:10.1016/S0169-5347(00)89061-8.21237012
Via S , LandeR. 1985. Genotype-environment interaction and the evolution of phenotypic plasticity. Evolution. 39 :505–522. doi:10.2307/2408649.28561964
Wagner A . 1994. Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization. Proc Natl Acad Sci USA. 91 :4387–4391. doi:10.1073/pnas.91.10.4387.8183919
Wagner A . 1996. Does evolutionary plasticity evolve? Evolution. 50 :1008–1023. doi:10.2307/2410642.28565284
Watson RA , WagnerGP, PavlicevM, WeinreichDM, MillsR. 2014. The evolution of phenotypic correlations and “developmental memory”. Evolution. 68 :1124–1138. doi:10.1111/evo.2014.68.issue-4.24351058
