On the evolutionary origin of discrete phenotypic plasticity

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.


Introduction
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.
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 , 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(Wagner , 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 E 1 and E 2 , 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 E 1 and inactive in E 2 .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.g i 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 (x i (0) = g i ).
x i changes through development, and let x i (τ) 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(τ) = (x 1 (τ), x 2 (τ), . . ., x n (τ)) during adult phase (i.e.τ v ≤ τ ≤ τ l ).The types-B and -C genes are involved in this developmental process.b ij 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 (c i 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: 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 x i 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 E 1 and δ(τ) = 0 in E 2 ).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 E 1 and E 2 at an interval of I T 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 E 1 if grown in E 1 and fit to E 2 when grown in E 2 from the same genotype (as illustrated in Supplementary Fig. S1).Let  X i = (X i1 , . . ., X in ) be the optimal expression levels of the n type-G genes in environment E i (X ij ∈ [0, 1/α]).The fitness of an adult with x(τ) in environment E i at developmental step τ is assumed to be given by exp Then, the lifetime fitness is defined as its geometric mean during the adult phase: 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 g i = 1/(2α), b ij = 0, and c i = 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 g i by Δg, where Δg is a random value drawn from a uniform distribution U( − γ g , γ g ).We assume 0 ≤ g i ≤ 1/α, so that Through a simulation run, we compute both w 1 and w 2 every generation, where w i is the fitness of the adult individual if it has grown in environment E i .The population averages are denoted by ̅ w 1 and ̅ w 2 , respectively.If the environment is E 1 , only w 1 is involved in selection and w 2 is evolutionary meaningless, and vice versa.Nevertheless, it is interesting to monitor both w 1 and w 2 through a simulation run to measure the potential ability of the system to fit the two environments.

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 b ij = 0, therefore the fitness (w 1 and w 2 ) is relatively low.When both w 1 and w 2 increase and become stably close to 1, we consider that phenotypic plasticity is attained.In each run, we simulated at most 10 7 generations and terminated when phenotypic plasticity successfully evolved.We considered that phenotypic plasticity had evolved at generation T = T p if the population mean fitness in both environments ( ̅ w 1 and ̅ w 2 ) was >0.95 in the following 10 5 generations.If phenotypic plasticity was not acquired until T = 10 7 , we considered that plasticity has not evolved successfully.In the following, we consider the parameter values listed in Table 1 unless otherwise specified.(0, 0, 5, 5, 0,  0, 5, 5) Origin of Phenotypic Plasticity | 3 Figure 2 summarizes the simulation results for various pairs of mutation rates and I T 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, S p , and T p 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 (T p ) tends to be short at intermediate I T (I T = 100, 1,000).It can be explained by the difference of the evolutionary dynamic between small and large I T .In the cases of a very small I T , the network should perform a challenging task to adapt to both environments almost simultaneously.For a large I T , the network can adapt to the environments in turn, perhaps making the adaptation process straightforward.However, there is a drawdown for a large I T .When I T 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 I T , such as I T ≥ 10 4 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.
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 E 1 and  X 2 = (0, 5) in E 2 .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 I T is very small and when I T is very large, therefore we consider two cases, I T = 1 and 10,000.The adaptation process to   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).
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 I T = 1, where the environment changes every generation and the system has to adapt to both environments almost simultaneously.Figure 3  ̅ 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).(x 1 (τ), x 2 (τ)) starts at × at embryo and moves along the red and blue arrows in E 1 and E 2 , respectively.The star represents (x 1 (τ c ), x 2 (τ c )) when the activation of the switch ends in E 1 .The default parameter values were used (see Table 1) except for ,  X 1 = (5, 0), and  X 2 = (0, 5).
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 b ij , c i , and g i 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 E 2 , (x 1 (τ), x 2 (τ)) moves directly towards the purple circle (as indicated by the blue arrow), while in E 1 , the evolved switch causes (x 1 (τ), x 2 (τ)) 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 (x 1 (τ l ), x 2 (τ 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 (x 1 (τ), x 2 (τ)) 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 E 2 , (x 1 (τ), x 2 (τ)) directly converges to the equilibrium near the optimal blue circle.Conversely, in E 1 , where the switch is activated, (x 1 (τ), x 2 (τ)) 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 (x 1 (0), x 2 (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 I T .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 E 1 at time point ▽ T1a, the equilibrium approaches (5, 0), indicating adaptation solely to E 1 .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 E 2 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 E 1 .Although the developmental outcome remains consistent across environments, we observe an additional equilibrium near the optimal expression for E 2 (i.e. the top-left corner).This new equilibrium is a consequence of the system's prior adaptation to E 2 , which persists even in the current phase within E 1 , 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 (x 1 (τ), x 2 (τ)) 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 E 1 and E 2 .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 (x 1 (τ), x 2 (τ)) 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 I T = 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.
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 b ij = 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 b ij , c i and g i to make adult phenotype close to its optimum in both environments.At T = 2 × 10 5 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 × 10 5 (see Fig. 5a).At time T = 2 × 10 6 (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 E 1 .In E 2 , 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 I T 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 I T .Figure 6 shows the result of a representative run when I T is very large (I T = 10,000).The pattern is quite different from that for a small I T 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 E 1 to E 2 (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 E 1 ), 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).
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 E 1 ), the system is fully adapted to E 1 , but not yet to E 2 .The expression patterns of the type-G genes in E 1 and E 2 are very similar to each other and The red and blue lines are the fitness in environment E 1 and that in E 2 , 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 × 10 5 , and 2 × 10 6 ) 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).
Origin of Phenotypic Plasticity | 7 fit to only E 1 .After I T = 10,000 generations at time point ▽ T1b (the end of E 2 ), the expression patterns of the type-G genes in E 1 and E 2 are again very similar to each other and fit to only E 2 .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 b ij , c i , and g i , 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 × 10 5 .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 b ij , c i and g i .The types-B and -C genes have contributed to this acquisition of plasticity, which is seen in well activated b ij and c i 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 I T where the system attempts to adapt to the two environments simultaneously, indicating I T 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 I T .We found that, in cases with intermediate I T , 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 S p and T p .In general, they have small effects on S p , while T p 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, T p decreases as μ c increases, irrespective of I T .The mutation rate of the type-B genes (μ b ) has a more complicated effect: μ b decreases T p when I T is small, while it slightly increases T p when I T is large.This pattern is explained by the difference in the evolutionary dynamics between the small and the large I T cases (see Figs. 5 and 6).When I T 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 T p .In contrast, when I T 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 Temporal change of fitness through evolution.The red and blue lines are the fitness in environment E 1 and that in E 2 , 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 × 10 5 , 3 × 10 5 , and 2 × 10 6 ) 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).
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 I T is large, μ b has a slightly positive correlation with T p .The effect of μ g seems small, suggesting that g i 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 T p and also have some impact on S p .As γ c increases, plasticity is more likely to evolve in a short time.A large γ b tends to decrease T p when I T is small, while increases T p when I T is large.S p becomes small in the parameter space where T p is very long (i.e.small γ c , and large γ b in large I T ).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, E 1 and E 2 , 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 E 1 and inactive in E 2 , or the expression level is 1/α and 0 in E 1 and E 2 , 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 E 1 and E 2 .We here consider intermediate environments between E 1 and E 2 , 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 E 1 to E 2 ), 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 w 1 and w 2 : When the distribution is discrete, in almost all range, either w 1 or w 2 is nearly 1, with only a small range where both w 1 are w 2 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, w 1 and w 2 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 E 1 or E 2 ).
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 w 1 and w 2 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 I T 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.

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 I T , as shown in Fig. 5), the network initially evolves to produce a phenotype with maximum "average" fitness (i.e. ������� w 1 w 2 √ ).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 coevolution 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 I T 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 I T 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 S p 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.
smaller than 0. The mutation rate μ b per generation applies to all b ij of the type-B genes, and each mutation changes b ij by Δb where Δb is drawn from U( − γ b , γ b ).The mutation rate of each c i 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.

Fig. 2 .
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 (S p ) and waiting time (T p ) are shown for various I T .In each panel, S p is shown in the upper bar plot while the distribution of T p 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 Table1).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).

Fig. 3 .
Fig. 3. Evolution of phenotypic plasticity in the simplified model with two type-G genes when I T = 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 E 1 and that in E 2 , 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).(x 1 (τ), x 2 (τ)) starts at × at embryo and moves along the red and blue arrows in E 1 and E 2 , respectively.The star represents (x 1 (τ c ), x 2 (τ c )) when the activation of the switch ends in E 1 .The default parameter values were used (see Table1) except for μ i = 10 −4 , σ = 5/

Fig. 4 .
Fig. 4. Evolution of phenotypic plasticity in the simplified model with two type-G genes when I T = 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 E 1 and that in E 2 , 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).(x 1 (τ), x 2 (τ)) starts at × at embryo, and moves along the red and blue arrows in E 1 and E 2 , respectively.The star represents (x 1 (τ c ), x 2 (τ c )) when the activation of the switch ends in E1.The default parameter values were used (see Table1) except for μ i = 10 −4 , σ = 5/

Fig. 5 .
Fig. 5. Evolution of phenotypic plasticity in the full model when I T = 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 E 1 and that in E 2 , 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 × 10 5 , and 2 × 10 6 ) 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 Table1).

Fig. 6 .
Fig.6.Evolution of phenotypic plasticity in the full model when I T = 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 E 1 and that in E 2 , 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 × 10 5 , 3 × 10 5 , and 2 × 10 6 ) 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 Table1).

Fig. 7 .
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 E 1 and E 2 b) are plotted.Default parameters and I T = 1 were assumed.c, d) An example case with a less discrete distribution.Default parameters except for μ c = 10 −6 and I T = 1 were assumed.

Table 1 .
Default parameter values in the simulation under the full model with eight type-G genes.
2 Optimal expression pattern in E 2