Evolution of Phenotypic Plasticity under Host-Parasite Interactions

Robustness and plasticity are essential features that allow biological systems to cope with complex and variable environments. Through the evolution of a given environment, the former, the insensitivity of phenotypes, is expected to increase, whereas the latter, the changeability of phenotypes, tends to diminish. However, in nature, plasticity is preserved to a certain degree. One possible cause for this is environmental variation, with one of the most important “environmental ” factors being inter-species interactions. As a first step toward investigating phenotypic plasticity in response to an ecological interaction, we present the study of a simple two-species system consisting of hosts and parasites. Hosts are expected to evolve to achieve a phenotype that optimizes fitness and increases the robustness of the corresponding phenotype by reducing phenotypic fluctuations. Conversely, plasticity evolves in order to avoid certain phenotypes being attacked by parasites. By simulating evolution using the host gene-expression dynamics model, we analyze the evolution of genotype-phenotype mapping. If the interaction is weak, the fittest phenotype of the host evolves to reduce phenotypic variances. In contrast, if a sufficient degree of interaction occurs, the phenotypic variances of hosts increase to escape parasite attacks. For the latter case, we found two strategies: if the noise in the stochastic gene expression is below a certain threshold, the phenotypic variance increases via genetic diversification, whereas above the threshold, it is increased due to noise-induced phenotypic plasticity. We examine how the increase in the phenotypic variances due to parasite interactions influences the growth rate of a single host, and observed a trade-off between the two. Our results help elucidate the roles played by noise and genetic mutations in the evolution of phenotypic plasticity and robustness in response to host-parasite interactions. Author summary Plasticity and phenotypic variability induced by internal or external perturbations are common features of biological systems. However, certain environmental conditions initiate evolution to increase fitness and, in such cases, phenotypic variability is not advantageous, as has been demonstrated by previous laboratory and computer experiments. As a possible origin for such plasticity, we investigated the role of host-parasite interactions, such as those between bacteria and phages. Different parasite types attack hosts of certain phenotypes. Through numerical simulations of the evolution of host genotype-phenotype mapping, we found that, if the interaction is sufficiently strong, hosts increase phenotypic plasticity by increasing phenotypic fluctuations. Depending on the degree of noise in gene expression dynamics, there are two distinct strategies for increasing the phenotypic variances: via stochasticity in gene expression or via genetic variances. The former strategy, which can work over a faster time scale, leads to a decline in fitness, whereas the latter reduces the robustness of the fitted state. Our results provide insights into how phenotypic variances are preserved and how hosts can escape being attacked by parasites whose genes mutate to adapt to changes in parasites. These two host strategies, which depend on internal and external conditions, can be verified experimentally, for example, via the transcriptome analysis of microorganisms.


Introduction
Robustness and plasticity are two important properties of biological systems. To 2 maintain function and high fitness in response to internal noise, environmental 3 variation, and genetic changes, the fitted state must be robust, whereas the phenotype 4 needs to possess plasticity in order to adapt to environmental variation. Indeed, the 5 evolution of robustness and plasticity has been investigated extensively both 6 theoretically and experimentally [1][2][3][4][5]. 7 In particular, the relationship between phenotypic fluctuations and robustness or 8 plasticity has recently attracted much attention [6]. As a result of evolution under a 9 given environment, phenotypic fluctuations generally decrease. With this decrease, the 10 system deviates less from the fitted state, thereby increasing robustness [7-9]. 11 However, decreased phenotypic fluctuation is also accompanied by decreased 12 adaptability to environmental variation, is discussed by generalising the 13 fluctuation-response relationship in statistical physics [6,10,11]. Here, phenotypes are 14 subject to the dynamic processes of several gene-dependent variables (for example, the 15 expressions of proteins). Thus, phenotypes vary according to genetic changes, 16 providing the phenotypic fluctuations [12][13][14]. On the other hand, as the dynamic 17 processes that shape phenotypes are influenced by external and internal noise, 18 phenotypes can vary even for isogenic individuals [15][16][17], representing another source 19 of phenotypic fluctuations [18,19]. Both these fluctuation sources depend on 20 genotype-phenotype mapping [20][21][22]. The evolution of genotype-phenotype mapping 21 is essential to understand the evolution of robustness and plasticity, and has been 22 explored extensively both in numerical [4,6] and laboratory experiments [23,24]. 23 Indeed, experiments and simulations of adaptive evolution under fixed conditions 24 have shown that fluctuations decrease over the course of evolution [25], and the 25 evolvability, that is, the rate of increase in the fitness or the phenotypic change per 26 generation, declines. Accordingly, as evolution progresses, the robustness of the 27 phenotype increases, while the plasticity decreases. 28 However, in nature, phenotypic fluctuations and evolutionary potentiality persist. 29 Even after evolution, the phenotype in question is not necessarily concentrated on its 30 optimal value, but its variance often remains rather large. For instance, even for 31 isogenic cells (clones), fluctuations in the concentration of each protein remains 32 sufficiently high to preserve potentially of evolution. This leads to ask the following 33 questions to be addressed: Why are such fluctuations not reduced, which, in principle, 34 would be possible by evolving appropriately negative feedback processes for 35 stabilization? How are phenotypic fluctuations or plasticity sustained? 36 One possible cause for the preservation of plasticity or fluctuation is environmental 37 fluctuation [26][27][28]. In natural evolution, fitness and environment are not fixed, but 38 variable. If environmental fluctuation is reduced excessively, the plasticity to adapt to 39 new conditions imposed by environmental changes would be lost. The plasticity of a 40 biological system dictates its capacity to cope with environmental changes that alter 41 fitness conditions. 42 For every species, interactions with other species are one of the most important 43 environmental factors that affect its existence. Even if a certain external 44 environmental condition itself is fixed, the types and populations of other species may 45 change owing to species-species interactions [29]. For instance, in the interaction 46 between hosts (prey) and parasites (predators), the former may change their 47 phenotype to escape attack by the latter, which, in turn, will change the phenotype of 48 the latter to continue attacking the former [30,31]. Thus, each species may retain 49 evolutionary plasticity to cope with dynamic changes in other species. Furthermore, 50 the phenotypic plasticity in one species may influence other species, resulting in a 51 dynamic change in inter-species interactions. If sufficient species-species interaction 52 occurs, the dynamic variation of phenotypes may be mutually sustained across 53 multiple species.

54
As a first step to investigate the phenotypic plasticity originating from ecological 55 interactions, we present a study of a simple two-species system consisting of hosts and 56 parasites. Within this system, distinct parasite types exist that can attack hosts 57 possessing a specific phenotype, whereas the host phenotypes are determined by 58 genotype-phenotype mapping designed to incorporate possible fluctuations. Hosts are 59 expected to evolve to achieve a phenotype that optimizes fitness and increases the 60 robustness of the successful phenotype by reducing phenotypic fluctuations. In 61 addition, they are also expected to evolve plasticity, that is, phenotypic adaptability to 62 cope with parasite attacks. We use a host gene-expression dynamics model to study 63 the evolution of genotype-phenotype mapping and the resultant phenotypic To study the evolution of genotype-phenotype mapping, we employed a simple 79 model for gene expression dynamics based on a sigmoidal function. This model 80 comprises a network consisting of genes that mutually activate or inhibit each 81 expression. In this model, there are M genes whose gene expression level, x (l) i 82 (i = 1, 2, ..., M ), corresponding to the host type l at time t is described by where J ij is the matrix representing the influence of gene j on the expression of gene i. 84 January 7, 2021 3/16 When j activates (represses) the expression i, it takes the form J ij = 1(−1), where 85 J ij = 0 if j does not influence i. The sigmoidal function, f (z), is given by for which β = 25, meaning that f (z) closely resembles a step function. Furthermore, 87 θ i represents the threshold of the input required for the expression of gene i. 88 Therefore, whether the gene is expressed (f (z) ≈ 1) or not (f (z) ≈ 0) depends on 89 whether the sum of inputs from the expression of other genes is larger than θ i .
which represents the stochasticity. The value of σ represents the strength of this noise. 99 Initially, all the gene expression levels are set smaller than the threshold θ i .

100
Furthermore, ϵ denotes the spontaneous expression level, which is smaller than θ i , and 101 dictates that x i must receive external or internal inputs from other x j values in order 102 to advance beyond θ i . I j is the external input for input genes, with I j = 1 for  Additionally, we set N = 300, M = 64, l inp = 8, and l out = 8.

106
Fitness and reproduction 107 The fitness of the host alone i.e., in the absence of parasites is determined by the 108 fraction of output genes that are expressed, as ds where x j is the expression of the output genes and µ 0 denotes the maximum fitness.

110
Therefore, the state with all output genes expressed (x j = 1, j = 1, 2, ..., l out ) is an 111 optimal phenotype for achieving maximal growth.  Below, we term the l p genes as target genes. Here, the target condition is given by 118 s(x i ) = 1 for x i > 0.5 and s(x i ) = 0 otherwise (as the threshold function is close to a 119 step function, x i takes ∼ 0 or ∼ 1 any way). For instance, a host whose expression 120 pattern (s(x 1 ), s(x 2 ), s(x 3 )) ≈ (0, 1, 0) is attacked by the parasite (i 1 , i 2 , i 3 ) = (0, 1, 0).

121
The growth rate, µ i , of each host decreases owing to the impact of the parasite i as 122 follows: where P m is the population density of the m-th parasite (m = 1, ..., N p = 2 lp ) and c is 124 the strength coefficient of the attack by the parasite. The volume of the host cells 125 grows according to When the volume exceeds a threshold value of 2 (v(0) set to 1), the host cell is 127 assumed to divide into two parts. If cell division leads to the cell number exceeds the 128 upper limit N by the cell division, the surplus cells are eliminated at random, in order 129 that the upper limit of the total cell number is N . Here, after division, the initial 130 expression of x i is reset to take a random value smaller than the threshold θ i .

131
Mutation is factored into the division process and added to the network J ij with a low 132 mutation rate. An (i, j) path in the network matrix is selected at random and its 133 values changed to one of the other to preserve the total path number. For example, if 134 J ij = 1(−1) and J kj = 0, then they are changed to J kj = 1(−1) and J kj = 0. We 135 prohibited direct connections between the "input" and "output" genes.
where κ ij takes unity only if the binary string of j can change to that of i via a single 148 mutation, otherwise it equals zero. To avoid the divergence of parasite populations, we 149 assume competition within all the parasite species, such that the total number of 150 populations is assumed to be fixed. Accordingly, we normalize the density P i such    Table 1).  222 We computed V ip and V g separately (see M ethods). The evolutionary time courses of 223 V ip and V g are shown in Fig 3. Here, the variances V ip and V g corresponding to all 224 output genes and V t ip and V t g of the target genes are plotted in Fig 3(ii). The cause phenotypic diversification (Fig 2). Figure 4 shows that the values of both V ip 241 and V g are low relative to the case with c > c t and σ < σ c . As V g is small, the 242 phenotype is not significantly changed by genetic variation. V ip and V g both decrease 243 while satisfying the inequality V ip > V g , leading to the evolution of robustness to noise 244 and genetic change. This supports previous studies [6,9,10]. The host population is of the other output genes are also maintained at high levels (see S1 Fig). Accordingly, 267 the fitness is reduced.

268
Transition between the Phase III and IV against noise strength 269 We also explored the transition between the last two phases. Figure 5 shows the Coping with parasites by phenotypic plasticity 279 In the absence of parasites, the variances V ip and V g decrease as the host 280 phenotype adapts to the environment, as described by the fitness condition through 281 evolution (Phase II). Next, we investigated whether introducing parasites after this 282 adaptation enables the evolution of plasticity to progress. Therefore, after the host 283 had adapted to the environment and acquired robustness (i.e., after V ip and V g had  . After this, the evolution of robustness was complete (V ip and V g decreases). Then, we introduced the interaction with the parasites by changing c from 0 to 3. The variances were computed from the isogenic variance over 20 iterations. The time course of (V t ip , V t g ) and (V ip , V g ) over generations after this switch is plotted for the target (circles; bold lines) and output genes (asterisks; dotted lines). The color bar indicates the number of generations since the switch. The variances V t ip and V ip both increase owing to interactions with the parasites, resulting in the evolution of phenotypic plasticity. Other model parameters: N = 300, M = 64, l inp = 8, l out = 8, and l p = 3.

Dependence of the variances on genetic change and by noise 293
In a previous study without host-parasite interactions [6], the correlated change 294 between V ip and V g during the evolutionary process was observed for σ > σ c . Under 295 the pressure of host-parasite interactions, the gene expression dynamics evolved to 296 diversify the phenotypes in the manner described above. Here, we discuss the increase 297 in V ip and V g through evolution. As shown in Fig 7(a), the two variance terms both variances maintain large values, so that the host can escape parasite attacks by 300 producing a different phenotype via mutation. In contrast , Fig 7(b) shows that the 301 variances of other output genes decreased and the evolution of robustness increased for 302 those phenotypes not attacked by parasites. However, this decrease is much smaller 303 than that corresponding to the absence of host-parasite interactions. To increase the 304 phenotypic variance of the target genes, the variances of other gene expressions need 305 to be maintained to a certain degree because of gene interactions through the GRN.

306
This trend of a weak increase in the phenotypic variance for all genes holds for c > c t 307 (see S2 Fig and S3 Fig).
for target genes (a) and the variance (V ip , V g ) for output genes (b) (c = 0 (asterisks) and 3 (circles)). The variances are computed from the isogenic variance over 100 iterations. The time course over generations is plotted for the variances of gene expression. The plots cover 3000 generations. We set the noise level to σ = 0.03 > σ c . In the absence of parasite interactions (c = 0), the two variances decrease, while V g < V ip is maintained throughout the evolutionary course. Conversely, under host-parasite interactions, V ip and V g show a correlated increase. Although this increase is prominent for target genes (asterisks in (a)), it is suppressed for output genes, before showing a slight decrease (b). 309 We examined how the increase in the phenotypic variances due to parasite interactions 310 influences the growth rate of a single host. Parasites concentrate their attacks on the 311 hosts with the most frequent phenotype in the population. Figure 8

319
In this study, we studied the evolution of phenotypic variances by using host gene 320 expression dynamics with the regulation network, in the presence of a host-parasite 321 interaction. If the interaction is weak, the host with the fittest phenotype evolves to 322 reduce phenotypic variances. In contrast, if the interaction is sufficiently strong, the 323 phenotypic variances of the host increase to escape specific phenotypes being attacked 324 by each parasite strain. We identified two strategies, either to increase the 325 noise-induced phenotypic variance (V ip ) or to increase the genetic variance (V g ), 326 depending on the strength of the noise in the stochastic gene expression. If the noise 327 strength is below the noise threshold, the diversification is primarily genetic in origin, 328 whereas above the threshold noise-induced phenotypic plasticity dominates, leading to 329 the variances V ip and V g increasing. In the latter case, both variances increase in 330 correlation, thus enhancing phenotypic plasticity, which helps to avoid parasite attacks. 331 Note that in a fixed environment without inter-species interactions, the two 332 variances V ip (due to noise) and V g (due to mutation) tend to decrease, thereby losing 333 phenotypic plasticity and evolvability. Under host-parasite interactions, the GRN 334 evolves to increase these variances, even after robustness has evolved and become  Whether biological populations deal with environmental changes by genetic 340 variation or phenotypic plasticity has been discussed both theoretically and 341 experimentally [3,32,33]. Indeed, the relevance of species-species interactions to 342 phenotypic plasticity and genetic diversification has been discussed extensively [34,35]. 343 Therefore, it is interesting to examine how the two strategies for phenotypic 344 diversification studied here are adopted therein.

345
First, we observed a trade-off between phenotypic plasticity and fitness. Hosts that 346 evolved plasticity to deal with parasite interactions tended to exhibit a decrease in 347 inherent fitness (i.e., the growth rate). Interestingly, such a trade-off has been 348 observed for predator-induced phenotypic plasticity [36][37][38] and in bacteria-phage 349 experiments [39], where the bacteria gain resistance to the phage by reducing the 350 competition for resources. Furthermore, it has also been suggested that excessively 351 strong phenotypic plasticity may reduce adaptability to climate change [40]. In 352 contrast, host populations that exhibit phenotypic diversification by increasing genetic 353 variation (V g ) in order to reduce the rate of infection do not show a significant decrease 354 in fitness. However, in this case, their phenotypes are less robust to noise or mutation. 355 Note that the generation of diverse phenotypes against uncertain environmental  Note that a stronger manifestation of genetic diversification results in speciation. In 381 our forthcoming paper, we aim to show how the phenotypic plasticity induced by 382 host-parasite interactions will also lead to genetic speciation.

384
We define the two phenotypic variances V ip and V g as follows.

389
(2) V (l) g (i) is defined as the phenotypic variance due to the genetic distribution.