Emergence of phenotypic plasticity through epigenetic mechanisms

Abstract Plasticity is found in all domains of life and is particularly relevant when populations experience variable environmental conditions. Traditionally, evolutionary models of plasticity are non-mechanistic: they typically view reactions norms as the target of selection, without considering the underlying genetics explicitly. Consequently, there have been difficulties in understanding the emergence of plasticity, and in explaining its limits and costs. In this paper, we offer a novel mechanistic approximation for the emergence and evolution of plasticity. We simulate random “epigenetic mutations” in the genotype–phenotype mapping, of the kind enabled by DNA-methylations/demethylations. The frequency of epigenetic mutations at loci affecting the phenotype is sensitive to organism stress (trait–environment mismatch), but is also genetically determined and evolvable. Thus, the “random motion” of epigenetic markers enables developmental learning-like behaviors that can improve adaptation within the limits imposed by the genotypes. However, with random motion being “goal-less,” this mechanism is also vulnerable to developmental noise leading to maladaptation. Our individual-based simulations show that epigenetic mutations can hide alleles that are temporarily unfavorable, thus enabling cryptic genetic variation. These alleles can be advantageous at later times, under regimes of environmental change, in spite of the accumulation of genetic loads. Simulations also demonstrate that plasticity is favored by natural selection in constant environments, but more under periodic environmental change. Plasticity also evolves under directional environmental change as long as the pace of change is not too fast and costs are low.


Introduction
Phenotypic plasticity is a phenomenon found in all kingdoms of life.It is defined as the ability of an organism to adjust its phenotype in response to the local environment and its changes, without any change in the genome (Pigliucci, 2005).Theoretical models suggest that phenotypic plasticity evolves under variable slow down evolutionary adaptation by genetic changes to permanently altered environmental conditions (Huey et al., 2003;Thibert-Plante & Hendry, 2011).Understanding fitness costs and limits of phenotypic plasticity requires considering the regulatory processes underlying phenotypic response patterns.Though some models cited above involve genetic-explicit individual-based simulations, the treatment of plasticity remains phenomenological, with no clear link to the action of the underlying plasticity mechanism.In nature, plasticity is typically mediated by regulatory processes (e.g., developmental switches, gene regulation, and neural or physiological control mechanisms, Sommer 2020).A change in plasticity, therefore, reflects a change in these regulatory mechanisms.
This change in perspective (from the phenomenological to the mechanistic approach to plasticity) has major implications (van Gestel & Weissing, 2016;Wagner, 2005;Watson & Szathmáry, 2016).The mechanistic approach has led to (i) many-to-one genotype-phenotype maps and, therefore, mutational robustness of the wild-type phenotype (van Gestel & Weissing, 2016), due to different genotypes potentially displaying similar phenotypic responses, and (ii) enhanced evolvability (i.e., ability to respond to selection, Houle 1992) by the uncovering of cryptic genetic variation when organisms are exposed to atypical environmental conditions.Cryptic genetic variation refers to genetic variation that has little or no effect on phenotypic variation under "normal conditions," but that under atypical (stressful) conditions generates phenotypic variation (Paaby & Rockman, 2014).However, it is unclear whether cryptic genetic variation results from hidden developmental programs that evolved from past selection or subject to drift when silent, or whether it is a manifestation of limitations of the plasticity mechanism itself, or a combination of both (Laland et al., 2015;Parsons et al., 2020;Romero-Mujalli et al., 2021).
In this paper, we propose a novel mechanism for phenotypic plasticity based on epigenetics, of the kind that may be attained via DNA-methylation/demethylation (Duncan et al., 2022).In nature, patterns of DNA methylation have been associated with both environmental factors and phenotypic variation in populations experiencing different-and stressful-environmental conditions, and observed for a wide range of taxa (e.g., bacteria, Casadesús, 2016;yeast, Hattman et al., 1978;snails, Fallet et al., 2020;Smithson et al., 2020;Thorson et al., 2017Thorson et al., , 2019;;plants, Niederhuth et al., 2016;mammals, including humans, Li & Zhang, 2014).Furthermore, the role of methylation on gene expression (silencing/activation) has been extensively documented and is phylogenetically widespread (Adrian-Kalchhauser et al., 2020;Curradi et al., 2002;Greenberg & Bourc'his, 2019;Kribelbauer et al., 2020).Epigenetic markers (e.g., methylation of the C5 positions in cytosines) can be induced by the environment and have been proposed as a candidate molecular mechanism underlying plastic phenotypic responses (Angers et al., 2020;Massicotte & Angers, 2012).However, the evolution and action of such a mechanism remains unknown since most of the empirical support is correlative (Richards et al., 2017;Schlichting & Wund, 2014).
In our model, the behavior (bind/release) of the epigenetic markers (methyl groups) follows random motion, but the frequency of events (methylation/de-methylation) is sensitive to organism stress (trait-environment mismatch).The epigenetic modifications (hereafter, epigenetic mutations or simply, epimutations) affect fitness-affecting polygenic traits of individuals.This is a novel approach that connects stimulus (stress) and response (trait expression) through feedback loops.On the one hand, this mechanism imitates developmental learning by trial and error, but on the other hand, it is also vulnerable to developmental noise leading to maladaptation.Using individual-based simulations, we study the role of epigenetic modifications on the performance of individuals and populations, as well as the evolution of the optimum level of plasticity under scenarios of periodic or directional environmental change.

Model and methods
We implemented an individual-based model where N individuals i = 1, … , N have phenotypes determined by a single response trait x i (t) that must match an environmental parameter p(t) across multiple generations t = 1, 2, 3, … .Natural selection favors the reproduction of phenotypes that are best at matching the environment.Response traits are genetically determined, but can be altered during a finite development period (there are two-time scales: developmental nested in generational).Such phenotypic plasticity (also called developmental plasticity; Botero et al., 2015) is enabled by epigenetic modifications, with adaptive rates  e controlled by a sensitivity trait  which also can be genetically coded, and thus evolve across generations.
Figure 1 illustrates the plasticity mechanism by the stress response of Theodoxus fluviatilis (Neritidae), a dioecious euryhaline snail that occurs in freshwater lakes and streams in central Europe.For T. fluviatilis, as for other snails under osmotic stress, p may represent the osmolality of the environment, and x the internal osmolality achieved by the accumulation of free amino acids and urea.Thus, x = p corresponds to a favorable isotonic state; x < p corresponds to a hypotonic condition (internal milieu less concentrated than external medium) leading to dehydration; and x > p to a hypertonic condition (internal milieu more concentrated than external medium) leading to swelling.
Our model implements these modifications by simulating DNA-methylation/demethylation that activate/deactivate multiple loci randomly, with alleles that determine response traits additively.This mechanism could also work for a single response trait locus, with an activation degree (0%-100%) given by a variable number of methylated CpG sites.Both interpretations are conceptually and mathematically identical, we keep the loci view because it is in line with quantitative-genetics theory and models (Barton et al., 2017;Bürger & Lynch, 1995;Scheiner et al., 2020).Once matured, individual traits determine mating success, and parents are replaced by offspring.Alleles for x i (t) and  i (t) pass to the next generation (t + 1) along with random mutations, but epigenetic modifications are lost, that is, no epigenetic inheritance (Jablonka, 2017).Natural selection acting directly on x i (t) and indirectly on  i (t) values takes place in environments with periodic or directional trends p(t) across generations.

Response and sensitivity traits
This section deals with independent development of single individuals in a generation.Thus, individual and generational markers (i, t) are omitted in favor of notation economy.An individual's response trait x is encoded by L diallelic loci as follows with allelic values l jk ranging from -∞ to ∞, that is, there are potentially infinite alleles per locus.k = 1, 2 denotes alleles from different parents.There are B non-plastic loci encoding a baseline trait value, and L -B plastic loci contributing or not, depending on activation bits or tags a j ∈ {1, 0}.The number of non-plastic B Figure 1.Conceptualization of our model, taking phenotypic plasticity in the snail T. fluviatilis as an example.An increase in the osmolality "p" of the environment (e.g., hypersalinity) leads to an osmotic imbalance associated with dehydration (tissue shrinkage).Stress ("x -p") sensing ("") modulate rates ("e") of epigenetic modification via DNA methylation (e.g., cytosine-methylation, C* in the DNA structure)/demethylation of multiple loci encoding proteins involved in the synthesis of metabolites (e.g., ninhydrin-positive substances, NPS).Increasing amounts of NPS represent the organism's response "x" which adjusts the internal osmotic situation of the animal to the external conditions, thereby stabilizing tissue volume.The organism's reaction alters its phenotype, but not its genotype.The photograph of a live freshwater snail has been taken by Amanda Wiesenthal.The originals of the redrawn graphs appeared in Symanowski and Hildebrandt (2010).The image of the DNA structure has been adapted from Hildebrandt et al. (2021).
loci must be at least equal to 1; otherwise, if B = 0, that is, all L loci are plastic, a scenario in which all loci become inactive would invoke an indeterminate trait value.This genotype-to-phenotype mapping implies that multiple genotypes can produce a common phenotype, but also a specific genotype can produce multiple phenotypes.To illustrate the first case suppose that a j = 1, then an individual can have a phenotype x = 2 if: (a) both alleles equal to one at one locus, all others zero; (b) one allele equal to one at two different loci, all others zero; (c) all alleles equal to 1/L; (d) large variation in allele values (some positive and some negative), with an average value of 1/L.To illustrate the second case, consider L = 10, B = 1, and l jk = 1: the phenotype x could be any even number from 2 to 20 depending on the activation states of L-B = 9 loci.Plastic loci j = B + 1, … , L experience epigenetic mutations over a development period of  days.These consist of random deactivations (a j = 1 → 0) and reactivations (a j = 0 → 1).The probability of epigenetic mutation per (plastic) locus per day increases with the mismatch between the response trait and the environmental parameter p (Foster, 2007).The sensitivity of the response is modulated by a second trait, .Thus, the response can be plastic: the phenotype of individuals with  > 0 can be modified during development by random loci activations and deactivations.
The sensitivity trait  is a critical parameter of (2).If it is neither too small nor too large, the frequency of epigenetic mutations is expected to decrease rapidly with |x -p|, making phenotype optimization by trial and error efficient.On the other hand, if  is too large, epigenetic mutations can be too frequent even for small |x -p|, resulting in developmental noise and maladaptation.Thus, optimum  values, that is, optimum levels of plasticity, are expected to result from natural selection.To account for this, we also consider  a heritable trait, encoded by a single locus as follows with allelic values  k also in the -∞ to ∞ range.A scaling parameter s > 0 is used to calibrate the simulations.

Population dynamics and evolution
Evolution of response x and sensitivity  traits occurs by natural selection acting on reproductive fitness, in populations with discrete, non-overlapping generations.A population of N individuals (i = 1, … , N) is split into M = 1/2N "males" (rounded to lowest integer) and F = N -M "females."Each female mates once; its mate is drawn at random from the population of males (with replacement, implying that males can mate more than once).The number of offspring produced by the pair is drawn from a Poisson distribution with parameter that is, the sum of female and male fitness.Hence, the "fitness" of an individual corresponds to the individual's effect on its expected number of offspring; this effect may include both survival and fecundity.We assume that male and female fitness is given by The term R exp (1 -N/K) may be viewed as a baseline fitness, which is density dependent: R corresponds to the "intrinsic reproduction factor," and the "carrying capacity" K is the population size above which baseline fitness drops below R.
factor in front of R introduces fitness depression experienced due to stress, which is defined as where x i is the response trait value (1) at the end of the development phase and p is the environmental parameter value in the current generation.Stress is 0 if x i = p, increases monotonically with the |x i -p| difference, and is bounded above by 1 (asymptote).
The inflection parameter  determines stress tolerance: small (large)  corresponds to low (high) tolerance, that is, strong (weak) natural selection.Plastic individuals bear a fitness cost C per unit of sensitivity trait ( i ).
Each offspring inherits one allele from each parent at each locus.Allele values of the inherited haplotype are randomly selected from the corresponding parental locus (no linkage).Allelic values (l jk ,  k , "Response and sensitivity traits" section) mutate with rate  m per locus per generation.Mutational effects (Δl jk , Δ k ) are normally sampled with zero mean and standard deviation  m .The offspring entirely replaces the parental population in the next generation.Keep in mind that since the offspring is the sum of a finite number (the number of females) of Poisson variables, extinction (the sum is zero) is possible.
Fitness is density-dependent in view of Equation ( 4), as well as frequency-dependent, because the fitness of couples (matings) combines heterogeneous male and female fitnesses.

Environmental change
We assume that individuals experience a constant environment p through their lives, but changes occur between discrete generations t = 1, 2, … We consider two trends: periodic or directional (Botero et al., 2015;Bürger & Krall, 2004).Cycles of period T and amplitude A are modeled by sinusoidal functions and directional changes are modeled by the linear function where  is the rate of environmental change.The initial state p 0 is zero'ed at the start of the simulations.Both trends are deterministic (i.e., no form of environmental stochasticity is considered).Although previous studies have shown that the color of the stochastic noise affects the importance of phenotypic plasticity for populations under environmental change (e.g., Ashander et al., 2016), its effect on phenotypic plasticity is beyond the scope of this study.

Simulations
We carried out the following simulation protocols for sexual diploid populations 1. Effect of plasticity on development.Follows the development of a population in a single generation.Individuals have distinct response traits x i updated according to (Equation (1)), but the sensitivity parameter  is the same for all (with  = 0 as control case).This allows us to study the action and limitations of the plasticity mechanism.2. Effect of plasticity on populations.Follows the evolution of response traits under allelic mutation and recombination.Like the first protocol, all individuals have the same fixed .This allows us to study how epigenetic-driven plasticity steers evolution under periodical or directional regimes of environmental change.Since individuals share the same , plasticity costs are meaningless here, so we set them C = 0. 3. Evolution of plasticity.Like the second protocol, alleles for the sensitivity trait also mutate, and  gets updated each generation according to (Equation (3)).Simulations start with  = 0 for all individuals, and we follow its evolution under different plasticity costs (C ≥ 0).We consider constant environment scenarios, as well as periodical or directional changes.
Response trait alleles l jk were initialized by normal sampling with zero mean and variance  2 G /2L, where  G = 1 is the standard deviation of traits in the population before development (first protocol) and evolution (second and third protocols).Sensitivity  1.

Effect of plasticity on development
We used single-generation simulations to demonstrate how epigenetic mutations enable a form of developmental learning that, by trial and error, adjusts phenotypes to current environmental conditions.Allelic values (l jk ) were normally sampled with mean zero and unit variance.Thus, response traits (x i ) are distributed approximately normally with a mean of zero at the start of development (notice probability density functions (p.d.f.) on the left sides of panels of Figure 2).All individuals share a common sensitivity trait value ().Plastic populations ( > 0) have a fair chance of matching phenotypes with the environment at the end of development (x mode ≈ p), as long as sensitivity traits are moderate.This is shown in Figure 2B and C as reduction of phenotype variance (right p.d.f.narrower than left p.d.f.), compared with nonplastic populations in panel A ( = 0, flat lines).However, if the sensitivity trait is too large like in panel D, plastic phenotypes are less likely to match the environment (final p.d.f. as wide as starting p.d.f.).Given the randomness of epigenetic mutations, phenotype plasticity does not progress in the same fashion for all individuals.Those starting with phenotypes close to the environmental parameter are unlikely to change at all during development, whereas those starting far from it display erratic trends toward the optimum phenotype value x i = p, or may miss the target by a large margin.Interestingly, if the sensitivity trait is extremely high, many individuals display phenotype oscillations and cannot approach the optimum.This is because epigenetic mutations rates (Equation ( 2)) can get very close to one, causing alternating activation/de-activation at plastic loci.These behaviors are illustrated in Figure 2, where we singled out the same four individuals ("cross", "circle," "triangle," and "square") in all panels.
To summarize, the effect of epigenetic mutations according to our model is not uniform.Optimal phenotypes are more likely to develop for intermediate values of the sensitivity parameter.If  or |x i -p| are too small, little to no activation changes take place; but if  or |x i -p| are too large, the high frequency of epigenetic mutations makes phenotype-environment matching unlikely.In other words, loss of function by developmental noise is pleiotropic with plasticity (and this limits the evolution of plasticity, see "Evolution of plasticity" section).Intermediate  values balance these tendencies: they make epigenetic mutations frequent enough at the start of development, which is useful to explore the phenotype space, and less frequent toward the end of development, preventing loss of optimality (Supplementary Figure A.1 shows a gradual decline of  e during development).We call this "learning-like behavior." Supplementary results show similar developmental trends if initial trait distributions are not normal, but asymmetric or bimodal (Supplementary Figures A.2 and A.3, respectively).These kinds of distributions could arise as a result of mutation and selection in past generations.We found that increasing the number of plastic loci (i.e., decreasing B, Supplementary Figure A.4) and the development time (increasing , Supplementary Figure A.5), raise the chances that optimal responses are found by developmental learning.

Effect of plasticity on populations
We wanted to assess the effect of plasticity on population dynamics and evolution.This was accomplished by simulating plastic ( = 0.03) and non-plastic ( = 0) populations during 2,000 generations.Genotypes are generated by recombination of parental genotypes from the previous generation with allelic mutation ( m = 10 -4 ).Evolution takes place under periodic (Equation ( 6), with period T = 200) or directional (Equation ( 7)) regimes of environmental change.The amplitude of periodic change was set at A = 1, and the rate of directional change at  = 0.001.This means that the amount of directional change matches the amplitude of periodic change in the 1,000th generation (or the periodic min-max range in the 2,000th).To explore the limits of phenotypic plasticity by epigenetic mutations, we also considered more extreme environmental changes.We multiplied the amplitude of periodic environmental variation and the rate of directional change by 4 and repeated the simulations (right column of Figures 3 and 4, respectively).Each {plastic, non-plastic} × {periodic, directional} combination was replicated 30 times.Figures 3 and 4 show the outcome of these simulations.
Population performance was assessed every generation by averaging the stress metric (5) of all individuals We also quantified the genetic load of the populations  1)) for 100 individuals (the same in all panels).The environmental parameter is p = 2, and the sensitivity trait is  = 0, 0.03, 0.1, 1, from panels A-D.The development of four individuals ("cross," "circle," "triangle," and "square") is highlighted.Trait distributions at the start and end of development are displayed on the left and right sides of each panel, respectively.Parameters from Table 1.
where W o and W o max are, respectively, the average and the highest fitness in a population without epigenetic alterations (i.e., a j = 1 in Equation ( 1)).A large genetic load means that a population consists of relatively many sub-optimal genotypes and relatively few genotypes with high fitness value (regardless of W o max being actually adequate or not in the current generation).
The efficiency of the plasticity mechanism was assessed in each generation by the ratio of phenotype variance between start (o) and end () of development V  /V 0 < 1 indicates that epigenetics drives a large number of genotypes toward common trait values, that is, phenotypes are canalized to match the environmental parameter p in the current generation.Increase in V  /V 0 ratios indicate that trait distributions are becoming a more wider or multimodal as a result of epigenetics.V  /V 0 > 1 indicates that epigenetic mutations are far too frequent ( e too large), leading to extreme phenotype variability.
When the environment changes periodically plastic and non-plastic populations alternate between low and high stress (Figure 3, panels A and B).Plastic populations tend to display lower stress levels compared with non-plastic, and differences are wider when the environment fluctuations are extreme.Indeed, under extreme periodic fluctuations, non-plastic populations experience stress levels that could be large enough to drive them toward extinction (panel B).Plastic populations accumulate larger genetic loads compared with non-plastic (panel C).This indicates that epigenetic mutations prevent elimination of large numbers of alleles by natural selection.This is possible because the environment changes between finite lower and upper bounds, and there is enough genetic variability in the population for plasticity to work upon.Indeed, when changes in the environment are moderate, phenotype variance ratios of plastic populations remain constrained over many generations (V t /V 0 < 1 on average, panel E), that is, populations are effectively canalized towards optimal phenotypes during periods of high and low stress.When changes are extreme though, genetic load in plastic populations attain very high levels (panel D), and the ability of epigenetic mutations to steer populations toward optimal phenotypes is lost after few generations (panel F).
Under directional change plastic and non-plastic populations maintain very low average stress levels (Figure 4, panels A and B).Speeding up environmental change (raising ) raises average stress levels for both population types.Genetic loads are low for both populations types but tend to be lower and decreasing for non-plastic (panels C and D).This indicates natural selection wiping out unfavorable alleles.Contrast this with plastic populations, where alleles can escape selection, allowing higher genetic loads (panel C), and its increase under rapid change (panel D).Like under periodic variation, plastic populations maintain low phenotype variance ratios, but there is an increasing trend nevertheless.This trend is weak under moderate environmental changes but becomes more evident when the rate of change is faster, leading to a sudden release of cryptic genetic variation, most of which is maladaptive (panel F), that is, trait adjusting by epigenetic mutations becomes less effective each generation.This likely happens because directional environmental change is asymmetric and unbounded: plastic populations with large genetic loads  8)), the middle row displays genetic loads (Equation ( 9)), and the bottom row shows phenotype variance ratios (Equation ( 10)).Ticker lines track averages.Parameters from Table 1.experience extreme environments with only a small fraction of individuals benefiting from epigenetic changes (those starting development close enough to optimum trait values), while the rest risk adaptation loss instead (due to developmental noise).

Evolution of plasticity
We simulated population dynamics with genetically encoded  (Equation ( 3)) in order to follow the evolution of plasticity itself.Without plasticity costs, sensitivity (plasticity) traits () evolve and tend to stabilize above zero in the long term.This happens in constant environments, for (moderate) periodic changes, and for (moderate) directional change (shown in Figure 5 by the average trend, for C = 0. Notice that the  attained after 2,000 generations is of the same order of magnitude as in Figure 2B and C).Further evolution of  halts because the benefits brought by epigenetic mutations, that is, developmental learning, are counteracted by the loss of function caused by developmental noise.
Plasticity costs are important in setting the extent of evolution: the higher the costs, the lower the long term  (Figure 6, see also temporal trends in Figure 5).Costs have much higher impacts under directional environmental change in comparison with periodic change.If the environment changes at a moderate pace,  values decline (near) uniformly with cost, but the drop is stronger under directional change compared with periodic (c.f., panels B and D).If directional changes are extreme, plasticity is not favored by selection, even at zero cost (panel E, compared with the neutral marker distribution).In stark contrast, under extreme periodic fluctuations sensitivity traits are favorably selected even at very high costs (panel C,  declines eventually with C, but is not shown).
We also simulated the evolution of  for asexual haploid populations (using simple versions of Equations (1 and 3) and without allele recombination).We found that the sensitivity (plasticity) trait is severely selected against in constant or directionally  8)), the middle row displays genetic loads (Equation ( 9)), and the bottom row shows phenotype variance ratios (Equation ( 10)).Ticker lines track average tendencies.Parameters from Table 1.

Discussion
This study demonstrates a novel mechanism in which epigenetic mutations enable phenotypic plasticity.This is achieved through random changes in the expression of multiple loci (or alternatively CpG sites in a locus), in a process that simulates the silencing/activation of genes caused by DNA-methylation (and demethylation).The frequency of epigenetic mutations correlates with the mismatch between the individual's phenotype and the external environment.That is, there is feedback between genetic expression and individual stress.As a result, a large set of heterogeneous genotypes can converge toward common phenotypes that better fit the environment.Epigenetic changes also depend on an individual's genotype; thus, plasticity itself can be subject to natural selection.
Our simulations predict that epigenetic mutations improve the match between individuals and the environment.This is also the case for populations under periodic or directional regimes of environmental change.However, the effectiveness of epigenetic mutations decreases if environmental changes are too extreme.This is because epigenetic alterations also promote large genetic loads (i.e., preservation of less favorable alleles), combined with the fact that our proposed mechanism is inherently imperfect: random epigenetic mutations cannot distinguish good from bad loci activations.This has important consequences for the evolution of plasticity.We found that even in the absence of explicit costs associated with the epigenetic mechanism, the evolution of plastic genotypes remains bounded: the benefits of phenotype space exploration by trial and error ("learning-like behavior") are balanced by the potential loss of function (developmental noise).1.

Learning-like plasticity
In our model, individuals, that is, genotypes, adjust their phenotypes to match an environmental parameter p (e.g., osmotic pressure) by experiencing it before natural selection takes place, that is, during a finite development time in which the random motion governing the behavior of epigenetic markers (implied to be methyl groups that bind/unbind parts of the genome) produces adaptive plasticity.In reality, binds (deactivations) and unbinds (activations) are carried out by methyltransferases and demethylases, respectively, enzymes capable of adding or removing  1. methyl radicals at specific places along DNA strands, affecting genetic expression (Moore et al., 2013;Shi et al., 2021).Thus, one could say that  controls the synthesis of these enzymes, or the proportion that are active through allosteric regulation.
The adjustment of the phenotype is unpredictable, in the sense that individuals do not "know" which parts of the genotype, that is, which loci, must be activated and de-activated in order to increase fitness.Thus, adaptation simulates a process of "learning by trial and error" but with "memory loss."A kind of "smart behavior" emerges due to a negative feedback on the rate of epigenesis (Equation ( 2)) during development: as the phenotypeenvironment mismatch |x -p| decreases, so does the frequency of epigenetic mutation  e , and the chances of losing an adapted phenotype are minimized.
The fact that randomness enables the "learning-like" process described above suggests that phenotypic plasticity could evolve without the (hypothetically) complex and costly biochemistry required to activate and de-activate "only the right loci."Thus, plasticity costs (Murren et al., 2015) can be very low, and random epigenetic mutations seem a plausible way to enable plasticity in simple, non-neural, unicellular organisms and in early life.It may also produce great morphological and physiological diversification in the face of highly conserved core genetic, cell biological, and developmental processes (Kirschner & Gerhart, 1998) (not tested in this study).Yet, this study shows that, in contrast to other models of plasticity (e.g., Lande 2009Lande , 2014;;Reed et al. 2010;Scheiner & Holt, 2012), perfect plasticity under negligible costs and high environmental predictability does not evolve.
For simplicity, our model considers a single phenotypic trait x and a set of loci that affect this trait either purely genetically (loci 1 to B) or by a combination of genetic and epigenetic effects (the remaining loci).We have shown that the epigenetic mutation rate can evolve to a value where a purely random process of placing epigenetic marks can markedly shift the trait x toward the optimal value p.The question arises as to whether this result still holds in the (more realistic) situation that (i) epigenetic mutations affect other genes (e.g., genes with "housekeeping tasks") for which plasticity is more harmful than beneficial and/or (ii) that different traits need to match different environmental targets.Regarding (i), the results of our study would, obviously, be much less convincing if the epigenetic mutation rate indiscriminately affect the functioning of all genes, including those that should stably encode the same phenotype, irrespective of the state of the environment.However, there is ample evidence that epigenetic marks are not placed indiscriminately and that the placement of these marks is related to the nucleotide sequence of the DNA (e.g.Adrian-Kalchhauser et al., 2020;Richards, 2006;Richards et al., 2017).In view of this, one would expect that natural selection on the DNA sequence will make "housekeeping" genes less susceptible to the placement of epigenetic marks.Regarding (ii), one solution is based on the fact that there is a diversity of epigenetic marks (e.g., Adrian-Kalchhauser et al., 2020).Different marks may be related to different "dimensions" of plasticity, that is, selection on the DNA sequence could make the loci encoding for trait x susceptible for epigenetic mark X and unsusceptible for mark Y, while the opposite happens for the loci encoding for trait y.Alternatively, several stress factors (perhaps even a large number) could together create an overarching level of stress (which might be a weighted average of the deviation of trait x from its optimal value p, of trait y from its optimal value q, etc.), which then, in turn, would regulate the epigenetic mutation rate for all loci determining x, y, etc.At present, much of this is speculative but an interesting topic for future investigations.

Genetic limits and emerging fitness costs to plasticity
Our study demonstrates two types of limits to phenotypic plasticity: (a) a genetic limit to plasticity, and (b) an intrinsic limit of the plasticity mechanism itself.
In nature, genomes are of finite size and this already imposes genetic limits to the range of phenotypes an individual organism can attain by means of its plasticity mechanism.Let us say that X is the set of possible phenotypes.Assuming x opt ∈ X, the intrinsic limit consists of the likelihood that x opt will be attained and not lost during the development phase due to epigenetic mutations.This depends on the sensitivity  (see Figure 2) and the traitenvironment mismatch |x -p|.Note that this limitation does not emerge because the plasticity machinery is energetically costly, but rather because the mechanism itself is vulnerable to catastrophic error if  is too high, or under extreme stress.In this case, developmental noise limits the ability of a plastic individual to match the optimal phenotype (DeWitt et al., 1998).
In this study, alleles constitute the raw material that the mechanism uses to develop phenotypic responses.According to our model, the genetic limits of plasticity can be expanded-though not necessarily always-by increasing the number of plastic loci L -B, for example by the insertion of transposable elements (Slotkin & Martienssen, 2007).In the model, each individual [genotype] has 2 L-B potential phenotypes thanks to epigenetic mutations (i.e., L -B binary states, active/inactive).Varying the number of plastic loci affected the degree of plasticity and effectiveness of the plasticity mechanism (see Supplementary Figures A.4,A.6,and A.7).Further study is needed to reveal why our model leads to a different conclusion than the study of Scheiner and Holt (2012) that found, that the number of non-plastic loci-but not the number of plastic ones-affects selection on plasticity.
Since the model uses real numbers for the allelic values, the condition all L loci active does not necessarily yields a maximum phenotypic value (e.g., maximum character size or length).However, it shows that the spectrum of phenotypic possibilities an individual can potentially attain during trait development, though diverse, is genetically limited by its genotype.It follows as well, that many loci of small effects can enable smoother phenotypic tracking of fluctuations in p, while few, stepwise, discrete responses.
Furthermore, it also demonstrates that different genotypes can have, to some extent, an intersecting set of phenotypic responses due to plasticity.The epigenetic mechanism can canalize such genetic variation into narrower phenotypic outcomes, which makes the optimal-wild type-phenotype mutationally robust (Hermisson & Wagner, 2004) (e.g., Figure 3C and E).The fact that genetically different individuals can still produce similar phenotypes enables the population to display coexistence of genetic polymorphisms, highly maintained genetic variation, and manyto-one genotype-phenotype maps (Wagner, 2005).Altogether, these features can enhance the evolvability of biological systems (Wagner, 2008).However, the coexistence of genetic polymorphisms due to the accumulation of mutations goes parallel with the accumulation of genetic load, and eventually with the loss of efficiency for the epigenetic mutation mechanism under strong regimes of environmental change (right column of Figures 3 and 4).Finally, the mechanism can make the same genotype susceptible to yield divergent phenotypic outcomes (e.g., when either,  or trait-environment difference |x -p| is too large, Equation ( 2)).
The model shows that plastic populations, when not able to minimize their mismatch with the environment (i.e., large |x -p|), release phenotypic variation, most of which is maladaptive.We observed this for directional trends, but may occur under periodic trends as well if the amplitude of fluctuations goes far beyond the population capability of producing suitable phenotypic adjustments (see Figures 3 and 4).Although the plasticity mechanism can effectively canalize populations toward optimal phenotypes, it inevitably leads as well toward increased genetic loads (see Figures 3 and 4).This variation is hidden, in the sense that it is already present in the population, and unleashed explosively (Figures 3F and 4F) due to large epigenetic mutation rates ( e ) failing to attain or keep suitable phenotypes.Under such extreme trends, plasticity hampers evolution.This provides an additional explanation to the often empirically observed cryptic genetic variation (Paaby & Rockman, 2014): a phenotypic outcome of the population (with high genetic load) that results from frequent failures of the plasticity mechanism in attaining the adaptive phenotype.Though to our model most of this variation is maladaptive, under other plasticity mechanisms (e.g., gene regulatory networks; van Gestel & Weissing, 2016), cryptic variation can be interpreted as opportunities potentially enhancing evolvability.

Evolution of the plasticity mechanism
The evolutionary simulations show that the plasticity mechanism introduced in this study can evolve in populations, starting from the non-plastic condition ( = 0).The evolution of plasticity (sensitivity trait) is bounded by developmental noise and by the cost of the plasticity mechanism.
Our simulations indicate that periodically fluctuating environments turn out to be the most favorable for the evolution of plasticity.This outcome may explain, for example, why the brackish water population of the snail T. fluviatilis displays a higher degree of plasticity than the freshwater counterpart.The salinity of the lakes inhabited by the fresh water populations is more constant along the year, as compared to the higher amplitude fluctuations of the salinity in brackish environments.
On the other hand and in agreement with previous theoretical studies, plasticity evolved under directional environmental change (Nunney, 2016;Scheiner et al., 2020).However, in contrast to our findings, under periodic environmental fluctuations, its evolution under directional environmental change was strongly impacted by plasticity costs and extreme environmental conditions.Thus, under the proposed mechanism, the potential importance of phenotypic plasticity in enhancing survival and in "buying time" for populations subject to directional-climate-change is questionable and demands empirical verification.
We performed our simulations under the assumption of sexual individuals, where genotypes are re-shuffled each generation due to allele recombination.Yet, our evolutionary simulations for haploids (see Supplementary Figure A.8) lead to the same general conclusion that plasticity is considerably more favored by natural selection under regimes of periodic rather than directional environmental change (or not changing at all).This could make our predictions testable by employing haploid laboratory models, for example, bacterial strains.

Conclusions
This study demonstrates that a classical explicit genetic model of quantitative genetics (Bürger & Lynch, 1995;Lynch et al., 1998;Scheiner et al., 2020) accounting for epigenetic mechanisms can produce adaptive, learning-like phenotypic plasticity; and, in-line with other mechanistic models (e.g., Brun-Usan et al., 2020a;van Gestel & Weissing, 2016;Wagner, 1994), developmental canalization of the genetic variation, many-to-one genotype-phenotype maps, mutational robustness of the wild-type phenotype, and the uncovering of cryptic variation.For previous mechanistic models, such attributes of the plasticity mechanism occur when assuming topological complexity in the form of number and types of regulatory connections (e.g., gene-regulatory networks; Brun-Usan et al. 2020b;van Gestel & Weissing 2016).This study expands this view and shows that they can as well result from random epigenetic mutations occurring on a trivial genotype-phenotype mapping (i.e., loci with additive interactions only).In addition, our mechanism produces: still functional phenotypic responses under novel environmental conditions (not previously screened by selection), genetic limits of plasticity, and emerging fitness costs to plasticity.All these features appear in the model under relatively simple underlying assumptions (i.e., random motion governing the regulatory process).These simple rules enable improving adaptation by a learning-like process.
Our mechanism is motivated by the following conceptualization of phenotypic plasticity: plasticity consists of the phenotypic change by which an organism attempts to return to a referential state (e.g., osmotic balance).Conserving such a state prevents function loss and, ultimately, favors increased fitness for plastic organisms.This is analogous to Le Chatelier's principle in chemistry, where "if a change is made to a system, the system will respond such as to absorb the force causing the change" (Anderson, 2005).
Given that the mechanism in this study requires developmental time and "trial and error," this study speculates that deterministic and exploratory mechanisms might be evolutionarily linked: faster and reliable deterministic mechanisms can accommodate (into complex topology) solutions found by exploratory mechanisms, which are typically prone to errors and may require longer developmental time.

Figure 2 .
Figure 2. Development of response traits x (Equation (1)) for 100 individuals (the same in all panels).The environmental parameter is p = 2, and the sensitivity trait is  = 0, 0.03, 0.1, 1, from panels A-D.The development of four individuals ("cross," "circle," "triangle," and "square") is highlighted.Trait distributions at the start and end of development are displayed on the left and right sides of each panel, respectively.Parameters from Table1.

Figure 3 .
Figure 3. Evolution of 30 plastic ( = 0.03) and 30 non-plastic ( = 0) populations under moderate periodic (left column, amplitude and period on top) or extreme periodic (right column, amplitude and period on top) environmental change.The top row displays population stresses (Equation (8)), the middle row displays genetic loads (Equation (9)), and the bottom row shows phenotype variance ratios (Equation (10)).Ticker lines track averages.Parameters from Table1.
changing environments, and favored if the environment changes periodically, in a similar fashion like Figure 6 (see Supplementary Figure A.8).

Figure 6 .
Figure 6.Sensitivity traits  attained by evolution for different cost levels.Boxplots show distribution of  averages for 100 populations after 2,000 generations.The distribution of a neutral marker (gray boxplot) serves as control.The environmental parameter p is constant in A; periodically changing (T = 200) in B (A = 1), C (A = 4); or directional in D ( = 0.001), E ( = 0.004).Parameters from Table1.

Table 1 .
Model variables and parameters.
k = 0 for the third protocol.The third protocol also considers a neutral marker (locus) with alleles initialized  k = 0, allowing us to distinguish the effects of natural selection from genetic drift.The allelic mutation rate ( m ) and effect ( m ) are the same for all allele types.A list of the model's variables parameters is shown in Table