Plasticity-led and mutation-led evolutions are different modes of the same developmental gene regulatory network

Environment

In plasticity-led evolution, the environment plays two roles: selector and inducer of phenotypes (West-Eberhard, 2003). We assume that these two roles of the environment are highly correlated (Ng & Kinjo, 2022; Ng & Kinjo, 2023).

We model the environment-as-selector as a 200-dimensional vector e of ±1. This vector is also the optimal phenotype used to calculate the fitness of each individual’s phenotype (also modeled as a 200-dimensional vector; see the following subsection). We model the environment-as-inducer as each individual’s environmental cue e. This vector is obtained by randomly flipping 5% of the elements of e.

We obtain a novel environment from a given ancestral environment by randomly flipping a given percentage of the elements of the ancestral environmental vector e. Note that a 50% change in the environmental vector is the maximum change that can be meaningfully achieved in our models. Each element takes +1 or −1 with equal probability. If we flip 50% of the elements of the environmental vector, then the correlation coefficient between the vectors before and after the change is 0. If we make a change greater than 50%, say 90%, then the correlation coefficient is −0.8, which is more significant than the case of 50% change. This implies that the GRN can flip all its state vectors (see below) and find a 90% match to the novel environment without introducing new mutations. In other words, GRNs are more adapted to a novel environment that is 90% different from the ancestral environment than one that is 50% different.

Development

In the following, we compare two variants (Ng & Kinjo, 2023) of the gene regulatory network (GRN) model introduced by Wagner (1996). In these models, the phenotype is a vector, and the genome is a set of interaction matrices. The first model, the Full model, incorporating environmental cues, developmental processes, and hierarchical regulation, exhibits PLE under large environmental changes (Ng & Kinjo, 2023). The environmental cue serves as the environmental input to the network. The developmental process is a recursive process that updates the system’s internal state, integrating environmental and genetic information into the phenotype. The hierarchical regulation reflects the multilayered nature of biological regulation.

In our developmental GRN models, let p(s) be the phenotype expressed at the s-th stage of development and $\tilde{p}\left(s\right)$ (and v(s)) be the exponential moving average (and variance) of the phenotype. To reflect the hierarchical regulation of GRN elements (such as epigenetic marks, RNAs, and proteins), we introduce layers of 200-dimensional vectors f(s), g(s),  and h(s) to represent epigenetic marks, gene expression and higher-order complexes (such as proteins, supramolecular complexes, etc.), respectively. For the Full model, we assume the following mutually recursive equations: (1)${\displaystyle f_i\left(s\right)={\sigma }_f\left(\sum _{j=1}^{200}{G}_{ij}{g}_j\left(s-1\right)+\sum _{j=1}^{200}{E}_{ij}\left(e_j-{\tilde{p}}_j\left(s-1\right)\right)\right),g_i\left(s\right)={\sigma }_g\left(\sum _{j=1}^{200}{F}_{ij}{f}_j\left(s\right)\right),h_i\left(s\right)={\sigma }_h\left(\sum _{j=1}^{200}{H}_{ij}{g}_j\left(s\right)+\sum _{j=1}^{200}{J}_{ij}{h}_j\left(s-1\right)\right),p_i\left(s\right)={\sigma }_p\left(\sum _{j=1}^{200}{P}_{ij}{h}_j\left(s\right)\right),{\tilde{p}}_i\left(s\right)=\alpha p_i\left(s\right)+\left(1-\alpha \right){\tilde{p}}_i\left(s-1\right),v_i\left(s\right)=\left(1-\alpha \right)\left\{v\left(s-1\right)+\alpha {\left[{\tilde{p}}_i\left(s-1\right)-p_i\left(s\right)\right]}^2\right\}}$ where E is a matrix that represents environmental regulation of epigenetic marks, F is a matrix representing epigenetic regulation of gene expression levels, G is a matrix representing genetic regulation of epigenetic marks, H is a matrix representing genetic regulation of higher-order complexes, J is a matrix representing interactions among higher-order complexes, P is a matrix representing regulation of the phenotype. The matrix ensemble E, F, G, H, J, P represents the individual’s genome. All matrices are sparse with a density of 0.02. The activation functions σ_f, σ_g, σ_h are based on the inverse tangent functions, and σ_p is based on the hyperbolic tangent function. These functions are modified in the spirit of LeCun (1989). The initial conditions are f(0) = 0, g(0) = 1, h(0) = 0 and p(0) = 0 where 0 is the zero vector and 1 is a vector where all elements are equal to 1. We recursively compute the state vectors f(s), g(s), h(s), p(s) for s = 1, 2, … until the total phenotypic exponential moving variance ∑_iv_i(s) < 10⁻⁵ for some s ≤ 200 (i.e., the phenotype has converged) and take ${\tilde{p}}_i\left(s\right)$ as the adult phenotype. Otherwise, we say that the phenotype does not converge.

The other model, the NoCue model, is identical to the Full model, except it uses the following equation for updating f_i(s) instead: (2)${\displaystyle f_i\left(s\right)={\sigma }_f\left(\sum _{j=1}^{200}{G}_{ij}{g}_j\left(s-1\right)-\sum _{j=1}^{200}{E}_{ij}{\tilde{p}}_j\left(s-1\right)\right).}$

The NoCue model does not take in environmental cues (e_j). As a result, it does not exhibit phenotypic plasticity and, therefore, serves as a model for MLE.

These models are mathematically equivalent to artificial recurrent neural networks, so they have generic learning capabilities. In this context, adaptive evolution can be viewed as a process of learning the environment by changing the genome (network weights) (Watson & Szathmáry, 2016). Notably, the Full model can adaptively respond to environmental cues through development.

Modeling evolution

We followed the same procedure as we did in Ng & Kinjo (2023) for simulating evolution by applying a genetic algorithm (Mitchell, 1998) to a population of 1000 individuals. We impose the role of the environment-as-selector by evaluating the fitness of each individual by matching the macro-environment e with the adult phenotype $\tilde{p}$ (assuming that $\tilde{p}$ converges). In other words, the macro-environment e is the optimal phenotype. We also included the number of developmental steps up to convergence into the fitness calculation so that individuals with fewer developmental steps are favored. We define the raw fitness of the ith individual as: (3)${\displaystyle {\omega }_i=exp\left(-\left(\alpha \parallel \tilde{p}-\mathbf{e}{\parallel }_1+\beta N_{\text{step}}\right)\right)}$ where $\parallel \tilde{p}-\mathbf{e}{\parallel }_1$ is the absolute (L1) distance between the first 40 (out of 200) elements of the adult phenotype $\tilde{p}$ and the corresponding elements of the macro-environment e, N_step is the number of developmental steps until convergence, α = 20 and $\beta =\frac{1}{20}$. Only 40 traits (elements of the phenotype vector) are subject to selection, and other 160(=200 − 40) traits are allowed to evolve freely. In this article, we call the value of $\parallel \tilde{p}-\mathbf{e}{\parallel }_1$ mismatch. Individuals whose phenotype $\tilde{p}\left(s\right)$ does not converge before a pre-specified number (200) of steps are given a zero fitness. The population in which all individuals have completed the developmental process is called the adult population in the following. The normalized fitness of the ith individual of the adult population is given by (4)${\displaystyle {\Omega }_i=\frac{{\omega }_i}{\underset{j}{max}\left\{{\omega }_j\right\}}}$ where max_j{ω_j} is the maximum raw fitness of the adult population. As such, Ω_i represents a probability of being selected for reproduction.

To generate offsprings, individuals are randomly sampled from the adult population with probability Ω_i. Each sampled individual is returned to the population before another individual is sampled. Selected individuals are paired to become parents. For each pair of parents, two offsprings are produced by randomly shuffling the corresponding rows of the genome matrices with a probability of 50% between the parents. The offspring population is duplicated to make two statistically identical populations of 1000 individuals each. In both populations, each element of the genome matrices of an offspring is independently and randomly mutated with a probability of 0.5% in a way that preserves the density of the genome matrices (Ng & Kinjo, 2023). We then let the two populations develop: one offspring population in the “ancestral” environment and the other in the “novel” environment. Only the individuals in the “novel” environment are subject to selection. The population in the “ancestral” environment is used only for comparison and is discarded after measurement. Adaptive evolution of a population is therefore modeled by repeated cycles of development, selection, and reproduction under the “novel” environment.

We let the population evolve in a constant environment e for 200 generations. This duration is called an epoch. Before simulating for another epoch, we change the environment e by flipping a specified proportion of the elements of this vector. We iterate this procedure for many epochs. Before the production runs presented in the Results section, we ran 40 epochs of preparatory runs under 40 completely uncorrelated environments to equilibrate the randomly initialized population and to let the system learn how to respond to environmental changes (Ng & Kinjo, 2023). In the following, the novel and ancestral environments refer to the current and last adapted environments, respectively.

Visualizing evolutionary trajectories

To visualize the evolutionary trajectory of a population, we project the phenotype and genotype of individuals in the population at each generation onto a two-dimensional genotype-phenotype space (Ng & Kinjo, 2023). First, the phenotype axis is defined as $\frac{{\mathbf{e}}_n-{\mathbf{e}}_a}{{∥{\mathbf{e}}_n-{\mathbf{e}}_a∥}_2^2}$ where e_n is the first 40 elements of the novel macro-environment, e_a is the first 40 elements of the ancestral macro-environment and ${∥{\mathbf{e}}_n-{\mathbf{e}}_a∥}_2$ is the Euclidean (L2) distance between e_n and e_a. We consider only the first 40 out of 200 elements because they correspond to the traits subject to selection. We project the phenotype p of each individual as: (5)${\displaystyle \mathbf{p}=\left(\tilde{p}-{\mathbf{e}}_a\right)\cdot \frac{{\mathbf{e}}_n-{\mathbf{e}}_a}{{∥{\mathbf{e}}_n-{\mathbf{e}}_a∥}_2^2}.}$ This way, the projected phenotypic values p of 0 and 1 correspond to phenotypes perfectly adapted to the ancestral and novel environments, respectively.

Next, the genotype axis is defined as follows. Let G_ij be the vector made by flattening and concatenating all the genome matrices of the ith individual of the jth generation. Let ${\overline{\mathbf{G}}}_j=\frac{1}{N}{\sum }_{i=1}^N{\mathbf{G}}_{ij}$ be the population average genotype vector of the jth generation. The genotype axis is defined as $\frac{{\overline{\mathbf{G}}}_{200}-{\overline{\mathbf{G}}}_1}{{∥{\overline{\mathbf{G}}}_{200}-{\overline{\mathbf{G}}}_1∥}_2^2}$. We project the genotype of the ith individual of the jth generation onto the genotype axis as: (6)${\displaystyle g_{ij}=\left({\mathbf{G}}_{ij}-{\overline{\mathbf{G}}}_1\right)\cdot \frac{{\overline{\mathbf{G}}}_{200}-{\overline{\mathbf{G}}}_1}{{∥{\overline{\mathbf{G}}}_{200}-{\overline{\mathbf{G}}}_1∥}_2^2}.}$ This way, projected genotypic values g_ij of 0 and 1 correspond to the average genotypes before and after one epoch of evolution in a novel environment, respectively.

We call the plot of the projected phenotypic values against the projected genotypic values the genotype-phenotype plot (Ng & Kinjo, 2023). Adaptive evolution generally proceeds from the lower left to the upper right on the genotype-phenotype plot.

Singular value decomposition of cross-covariance matrices

We use cross-covariance matrices to study the correlation between selected phenotypes and environmental noise or genetic variation. We define the Pheno-Cue or Pheno-Geno cross-covariance matrices as (7)${\displaystyle C_{ij}=\frac{1}{N}\sum _{k=1}^N\left(p_{ik}-{\overline{p}}_i\right)\left(x_{jk}-{\overline{x}}_j\right)}$ where p_ik is the ith trait of the phenotype vector of the k-th individual, ${\overline{p}}_i$ is the population average of the ith trait, x_jk is the jth environmental factor (for Pheno-Cue) or the jth element of the vectorized genome (for Pheno-Geno) of the k-th individual and ${\overline{x}}_j$ is the corresponding population average value.

We apply singular value decomposition (Yanai, Takeuchi & Takane, 2011) to the cross-covariance matrix C to obtain orthonormal components as follows (Navarra & Simoncini, 2010): (8)${\displaystyle C=\sum _i{\sigma }_i{u}_i{v}_i^{\top }}$ where the superscript ⊤ indicates transpose. In Eq. (8), u_i and v_i are the ith left and right singular vectors, interpreted as the ith principal axis of phenotypic variation (the left singular vectors) in response to the corresponding principal axis (the right singular vectors) of environmental noises (or genetic variation); σ_i is the ith singular value, arranged in decreasing order, which represents the cross-covariance between the left and right principal components.

To quantify developmental bias, we used the proportion of the (squared) first singular value to the total cross-covariance: (9)${\displaystyle \frac{{\sigma }_1^2}{\sum _i{\sigma }_i^2}.}$

The alignment between the first principal axis of phenotypic variation and environmental change is measured by the magnitude of the normalized dot product: (10)${\displaystyle \frac{|u_1\cdot \left({\mathbf{e}}_n-{\mathbf{e}}_a\right)|}{\parallel {\mathbf{e}}_n-{\mathbf{e}}_a\parallel }}$ where e_n is the novel environment and e_a is the ancestral environment.

The Full model exhibits plasticity-led evolution under large environmental changes

This subsection briefly summarizes our previous results (Ng & Kinjo, 2023) to introduce fundamental notions describing PLE and MLE. We simulated evolution under large environmental changes for ten epochs, flipping 50% of the elements of the environmental vector at the end of each epoch. We visualized evolutionary trajectories using the genotype-phenotype plots (Figs. 1A, 1B; see also Video S1) (Ng & Kinjo, 2023). Phenotypic values of 0 and 1 correspond to phenotypes perfectly adapted to the ancestral and novel environments, respectively. Similarly, genotypic values of 0 and 1 correspond to the average genotypes adapted to the ancestral and novel environments. Generally, the evolution proceeds from the lower left corner to the upper right corner of the genotype-phenotype plot.

For the full model (Fig. 1A), phenotypic values are greater in the novel environment than in the ancestral environment on average, indicating an adaptive plastic response. This response is also called phenotypic accommodation (West-Eberhard, 2003; West-Eberhard, 2005). The standard deviation in the phenotypic values in the ancestral environment is minimal, indicating the robustness of adapted phenotypes. The standard deviation in the phenotypic values in the novel environment is much greater than that in the ancestral environment in the first generation, indicating an increase in phenotypic variability. Our previous work shows that this is attributed to the uncovering of cryptic genetic variation (Ng & Kinjo, 2023). The phenotypic and genotypic values in the novel environment increase somewhat quickly, indicating rapid genetic accommodation. These features of the genotype-phenotype plot assert that the Full model exhibits plasticity-led evolution under large environmental changes.

The NoCue model (Fig. 1B) exhibits neither adaptive plastic response nor uncovering of cryptic genetic variation (by design). The cluttered distribution of points in the lower left corner demonstrates that the early adaptation stage of these trajectories is prolonged. The sparse distribution of points around the center of the plot indicates a rapid change in genotype once adaptive mutations are found.

We next tracked the trajectories of mismatch (Fig. 1C). Recall that the mismatch measures the adaptedness using the distance between the phenotype and environment (see the subsection ‘Evolution’ in ‘Materials and Methods’). The mismatch of the Full model decreases immediately and rapidly. In contrast, the NoCue model exhibits a delay of around 20–30 generations before the mismatch decreases rapidly, a manifestation of gradualism.

Finally, we checked the trajectories of genetic variance (Fig. 1D). Here, the genetic variance is computed as the sum of the variance in each element of the genome matrices. For the Full model, the genetic variance decreases immediately and rapidly to a minimum before gradually increasing again. The initial rapid decrease corresponds to a purifying selection of uncovered cryptic variations, whereas the later increase indicates the accumulation of cryptic genetic variation. This behavior contrasts the NoCue model, for which the genetic variance gradually increases at the start of each epoch before decreasing to a minimum, after which the genetic variance gradually increases again towards the end of the epoch. The initial increase in the genetic variance coincides with the initial delay of mismatch (Fig. 1C), which is attributed to the time taken to find adaptive mutations in the novel environment. We identify these initial behaviors of the NoCue model as a key signature of MLE.

The genetic variance of the Full model is consistently less than that of the NoCue model (Fig. 1D), indicating that the Full model harbors less standing genetic variation than the NoCue model. Nevertheless, the Full model adapts much more quickly than the NoCue model. This is because most genetic variation remains latent in the NoCue model due to the developmental system’s robustness and the absence of environmental cue input.

The Full model exhibits mutation-led evolution under small environmental changes

We simulated evolution under small environmental changes for ten epochs, flipping 5% of the elements of the environmental vector at the end of each epoch. The evolutionary trajectories of the Full and NoCue models (Figs. 2A, 2B) are nearly identical under small environmental changes. That is, the Full model (Fig. 2A) exhibits neither adaptive plastic response nor uncovering of cryptic genetic variation. The NoCue model exhibits similar evolutionary trajectories under small (Fig. 2B) and large (Fig. 1B) environmental changes. Therefore, under small environmental changes, both the Full and NoCue models exhibit the characteristics of MLE. This observation is consolidated by the trajectories of mismatch and genetic variance (Figs. 2C, 2D).

The Full model undergoes a sharp transition from MLE to PLE

We have seen that the Full model exhibits PLE or MLE under large or small environmental changes. This observation suggests a transition from MLE to PLE as environmental change increases. To see how this transition occurs, we examined the mean projected phenotype (Fig. 3A) in the first generation after environmental changes of different magnitudes. The mean projected phenotypic value suddenly increases (indicating adaptive plastic response) at around 30% environmental change. The variance in projected phenotype (Fig. 3B) also suddenly increases (indicating the uncovering of cryptic genetic variation) around 30% environmental change.

To further analyze the phenotypic variation under various environmental changes, we examined the cross-covariance matrix between phenotype and environmental cue (Pheno-Cue cross-covariance) or genetic variation (Pheno-Geno cross-covariance). These cross-covariance matrices quantify how phenotypic variation is affected by environmental noise or genetic variation. In particular, we only consider the first singular value of the Pheno-Cue (Fig. 3C) and Pheno-Geno (Fig. 3D) cross-covariance matrices as it explains most phenotypic variation. We found that the first singular values also suddenly increase around 30% environmental change. This observation indicates that the phenotypes become susceptible to environmental noise or genetic variation, suggesting a breakdown of environmental and mutational robustness around this point.

Developmental bias precedes PLE

The cross-covariance matrices introduced in the previous subsection contain various information about phenotype distribution, particularly the extent of phenotypic diversity and the skew in phenotype distribution (Noble, Radersma & Uller, 2019). There, we studied the extent of phenotypic diversity measured by the first singular values. We now study the skew in phenotype distribution, also called the developmental bias (Uller et al., 2018). In the previous work, we observed large developmental biases aligned with large environmental changes, whereas such biases were absent in adapted environments (Ng & Kinjo, 2023).

To study the transition in developmental bias, we measured the proportion of the first singular value of the Pheno-Cue (Fig. 4A) and Pheno-Geno (Fig. 4B) cross-covariance matrices over environmental changes of different magnitudes. There is an increase in both proportions at around 25% environmental change. Compared to the transition from MLE to PLE (Fig. 3), this increase in bias is less abrupt and occurs at smaller environmental changes. Thus, the developmental bias is more susceptible to environmental changes than the extent of phenotypic diversity. Still, it does not manifest as adaptive plastic responses or increased phenotypic variance due to the system’s robustness when the environmental changes are small. When the bias becomes sufficiently large, the system’s robustness breaks down, as indicated by an abrupt increase in the cross-covariances (Figs. 3C, 3D).

We next examined the orientation of the developmental bias. Specifically, we measured the alignment between the environmental change and the first phenotype singular vector of the Pheno-Cue (Fig. 4C) or Pheno-Geno (Fig. 4D) cross-covariance matrices over different magnitudes of environmental changes. The alignment is small for less than 25% environmental changes. Above 25% environmental change, developmental bias is well-aligned with environmental change. Thus, the emergent developmental bias facilitates selection by magnifying the difference between the most and least fit individuals. Note that the transition in alignment also happens before the increase in phenotypic variance (see Figs. 3C, 3D). Although plastic responses are invisible at this size of environmental change, such developmental bias aligned with the environmental change is expected to facilitate rapid adaptation should a large environmental change occur (Uller et al., 2018; Noble, Radersma & Uller, 2019).