Bastien Polizzi

Vincent Calvez

Sylvain Charlat
email: sylvain.charlat@cnrs.fr

Etienne Rajon
email: etienne.rajon@univ-lyon1.fr

(Epi)mutation rates and the evolution of composite trait architectures

Mutation rates vary widely along genomes and across inheritance systems. This suggests that complex traits -resulting from the contributions of multiple determinants -might be composite in terms of the underlying mutation rates. Here we investigate through mathematical modeling whether such an heterogeneity may drive changes in a trait's architecture, especially in fluctuating environments where phenotypic instability can be beneficial. We first identify a convexity principle, related to the shape of the trait's fitness function, setting conditions under which composite architectures should be adaptive or, conversely and more commonly, should be selected against. Simulations reveal, however, that applying this principle to realistic evolving populations requires taking into account pervasive epistatic interactions that take place in the system. Indeed, the fate of a mutation affecting the architecture depends on the (epi)genetic background, itself depending upon the current architecture in the population.

We tackle this problem by borrowing the adaptive dynamics framework from evolutionary ecologywhere it is routinely used to deal with such resident/mutant dependencies -and find that the principle excluding composite architectures generally prevails. Yet, the predicted evolutionary trajectories will typically depend on the initial architecture, possibly resulting in historical contingencies. Finally, by relaxing the large population size assumption, we unexpectedly find that not only the strength of selection on a trait's architecture, but also its direction, depend on population size, revealing a new occurrence of the recently coined phenomenon of 'sign inversion'.

Introduction

Complex traits, those resulting from the combined contributions of multiple determinants, may in principle be subject to composite mutations rates. For example, genes contributing to the same trait may mutate at rates spanning orders of magnitude because their genomic environments differ, or simply their lengths [START_REF] Rando | Timescales of genetic and epigenetic inheritance[END_REF][START_REF] Hodgkinson | Variation in the mutation rate across mammalian genomes[END_REF][START_REF] Oman | How Sequence Context-Dependent Mutability Drives Mutation Rate Variation in the Genome[END_REF].

The possibility of heritable yet non-genetic contributions introduces even more pronounced heterogeneity in mutation rates sensu lato, that is, in the rates of random change, whatever the inheritance mechanism [START_REF] Rando | Timescales of genetic and epigenetic inheritance[END_REF][START_REF] Danchin | Avatars of information: towards an inclusive evolutionary synthesis[END_REF]. Indeed, non-genetic inheritance systems, as multiple as they may be, share the property of being less stable than DNA sequences [START_REF] Johannes | Assessing the impact of transgenerational epigenetic variation on complex traits[END_REF][START_REF] Rechavi | Transgenerational Inheritance of an Acquired Small RNA-Based Antiviral Response in C. elegans[END_REF][START_REF] Denkena | Region-level epimutation rates in Arabidopsis thaliana[END_REF][START_REF] Dodson | Heritable capture of heterochromatin dynamics in Saccharomyces cerevisiae[END_REF][START_REF] Graaf | Rate, spectrum, and evolutionary dynamics of spontaneous epimutations[END_REF]. Here we investigate whether such an heterogeneity in mutation rates may drive the evolution of a trait's architecture: just as mutation rates can be adaptively tuned by the degree of environmental stability [START_REF] Ishii | Evolutionarily stable mutation rate in a periodically changing environment[END_REF][START_REF] Johnson | Beneficial mutations, hitchhiking and the evolution of mutation rates in sexual populations[END_REF][START_REF] Andre | The evolution of mutation rate in finite asexual populations[END_REF], can trait architectures evolve toward optimally mixed contributions of determinants differing in their rates of mutations? Or, on the contrary, are non-composite mutation rates generally selected for, which should lead to simplified traits architectures, homogeneous in terms of the underlying mutation rates?

We address these questions through mathematical modeling, assuming that selection acts on a single trait in an environment oscillating between two states. We further assume that the relative contributions of two (or more) determinants to the trait, differing only in their respective mutation rates, are controlled by one evolving parameter, and we study its evolution. In our analysis, the mutation rates differ by orders of magnitude, such that they can be referred to as 'fast' and 'slow' in the case of two determinants. We start with the analytical treatment of an idealized situation where only the absolute growth rate of a mutant determines its fate, and find the architectures that maximize this fitness proxy.

We thereby identify a general mathematical principle setting boundaries to the conditions where complex architectures, and thus composite mutation rates, may be optimal. According to this principle, whatever the degree of environmental instability, non-composite mutation rates are favored as long as individuals of intermediate phenotypic value do not perform better across different environments than the average of the extreme phenotypes. In other words, if the fitness function of the trait is convex, favoring specialists over generalists, the fitness function of the trait architecture is also convex such that non-mixed contributions are selectively favored. However, we also show that, as long as the degree of environmental instability is not too high, non-composite mutation rates remain optimal even if the fitness function of the trait is concave.

Absolute growth rates, however, ignore a potentially important aspect of the evolutionary dynamics at play in this system: a mutation necessarily occurs in an (epi)genetic background whose state (i) may depend on the strategy in place, and (ii) may critically affect the mutant's fate. We address this issue by using the adaptive dynamics framework, initially built to deal with ecological interactions by accounting explicitly for context-dependencies in a mutant's invasion success [START_REF] Metz | How should we define 'fitness' for general ecological scenarios?[END_REF][START_REF] Geritz | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF]. This approach generally confirms the above outlined convexity argument, but also reveals that pervasive epistatic interactions render the system reluctant to change, because determinants contributing little to the trait tend to accumulate cryptic variation that becomes deleterious if the trait architecture is changed. This reveals that historical contingencies may prove important in the evolution of genetic architectures.

Finally, we use simulations to relax the assumption of a large population size and thus investigate the consequences of genetic drift on the evolution of a trait's architecture. In doing so, we reveal a counterintuitive pattern whereby reducing population size below a threshold changes not only the efficiency but also the direction of selection. We relate this finding to a recently described phenomenon, coined 'signinversion', that generally takes place whenever some source of variability produces selective pressures differing in their direction, average duration, and strength [START_REF] Raynes | Sign of selection on mutation rate modifiers depends on population size[END_REF][START_REF] Raynes | When good mutations go bad: how population size can change the direction of natural selection[END_REF]. In such situations, weak but stable selection in one direction is dominant in large populations, while strong but occasional selection in the other direction becomes the main player in small populations. In our case, sudden environmental changes strongly select for a large contribution of the fastest determinant whereas environment stasis selects for a reduced input of mutations and therefore a high contribution of the slowest determinant. Yet, the latter selection pressure is typically weak and long-lasting, thereby acting efficiently in large populations only, resulting in a trait's architecture that is contingent on population size.

Models and tools

Qualitative summary of the model

Trait architecture. The model considers I ≥ 2 determinants with distinct mutation rates that contribute to a phenotypic trait, varying continuously from 0 to 1. Each determinant is characterized by (i) its (epi)genotype, taking a binary value (X i = 0 or X i = 1), (ii) its mutation rate (µ i > 0) at which the (epi)genotype switches randomly from X i to 1 -X i , and (iii) its relative contribution to the trait (α i ∈ [0, 1]). The contributions sum to one, that is i α i = 1. The resulting trait for a given individual is the weighted sum of its determinant values: Φ X (α) = i α i X i (simply denoted as Φ if unambiguous).

In the following, the vector of contributions α = (α i ) 1≤i≤I of a single individual is referred to as its trait architecture. In the simulations, we assume that the values of α are subject to mutational changes, at the slowest rate min µ i . Mutated α values are drawn in a Gaussian law centered on the parental value with variance 0.1, and constrained to satisfy i α i = 1.

Notations in the case of two determinants (fast versus slow). The main part of our study is devoted to the case of two determinants, below referred to as 'fast' and 'slow', in reference to their contrasted mutation rates µ F and µ S (µ F ≫ µ S ). In this particular case, we can denote the respective contributions of the fast and slow inheritance systems as α F and α S = 1 -α F . The (epi)genotype values are then denoted X F and X S , both in {0, 1}.

Evolutionary dynamics. We use a birth and death model, continuous in time, so that generations are overlapping. The fecundity rate is assumed to be homogeneous. Up to changing units of time, it is set to 1. Selection acts on the rate of mortality, namely s|Φ -Φ * (t)| γ , where the strength of selection is controlled by parameter s, Φ * is the environment-dependent optimal phenotype, and the exponent γ ≥ 0 is the shape parameter of the fitness function (see Figure 1). We distinguish three cases on the basis of γ values: γ = 1 is our reference case, and corresponds to a linear relationship between trait and fitness; γ < 1 corresponds to a convex relationship resulting in a higher average fitness of specialists (Φ close to extremal values 0 or 1); in contrast γ > 1 results in a concave relationship favoring generalists (Φ close to 1 2 ).

We assume that population size is regulated by an additional density-dependent contribution to the mortality rate N (t) K , where N (t) is the instantaneous population size, and K denotes the carrying capacity.

Environmental variation. We assume that the environment fluctuates between two possible states A/B associated with different optimal phenotypes Φ * A = 0 and Φ * B = 1. We consider periodic changes of period T , with a symmetrical alternance between A and B. It is relevant to distinguish three regimes of environmental variation in relation with mutation rates. First, the environment may be changing very rarely, that is, even more slowly than the slowest determinant (T (min µ i ) ≫ 1). Second, the environment may be changing very often, that is much faster than the least stable determinant (T (max µ i ) ≪ 1). A third, intermediate regime, is one where the rate of environmental change is comparable to the range of mutation rates.

Previous occurrence of the model. The model is inspired by the study in [START_REF] Rajon | In)exhaustible Suppliers for Evolution? Epistatic Selection Tunes the Adaptive Potential of Nongenetic Inheritance[END_REF] who restricted their model to the case of two determinants (I = 2), with respective contributions α 1 = 1 -α and α 2 = α, and to γ = 1. In their interpretation, the number α ∈ [0, 1] measures the balance between genetic and epigenetic contributions to the trait value. The only difference lies in the mutation rate, being significantly larger for the epigenetic contribution. While not limiting our interest to epigenetic inheritance, we retain from this initial work the focus on the variability of the mutation rates µ i 's and its consequences on the evolution of α. Separation of mutation timescales remains indeed the main component of our mathematical analysis.

A Lyapunov exponent approach in infinite populations

We begin our analysis by considering a monomorphic population for the trait architecture α = (α i ) 1≤i≤I , taking its growth rate as a proxy for its fitness. To this end, we introduce the Lyapunov exponent λ(α), which is the overall growth rate of the population over a complete cycle of environmental fluctuations.

The rate λ(α) is measured in the long run, that is, once the respective frequencies of the (epi)genotypes (X i ) 1≤i≤I have reached their periodic equilibrium distribution f g (t), where g denotes the (epi)genotype, taking value in {0, 1} 2 I . Technically, this equilibrium distribution is characterized by a periodic eigenvalue problem (Floquet spectral theory, see Appendix A), with the eigenvalue λ(α) defined by the following formula:

The Lyapunov exponent approach described above ignores the influence that the state of the population in place (denoted the 'resident', with trait architecture α r ), may have on the population dynamics of a 'mutant' with a different architecture α m . The adaptive dynamics framework allows to take these interactions into account by calculating the mutant's growth rate at the time of invasion -also known as invasion fitness -assuming that its own influence on the resident's population dynamics can be neglected at this initial stage.

In our model, the way the resident may exert an influence on the mutant's fitness is through its impact on the equilibrium frequencies of the (epi)genotypes, hence providing a favorable or unfavorable context for a different genetic architecture. For example, in the two-determinants case, and assuming a stable environment, if the resident population happens to be monomorphic for α F,r = 1 and thus α S,r = 0, then X S reaches its neutral equilibrium distribution (1/2; 1/2). This means that a mutant with non zero α S,m will on average have lower fitness than the resident, even though the absolute growth rate, as measured by the Lyapunov exponent approach, would be highest for α S = 1.

To take such resident / mutant interdependence into account, we computed in the two-determinants case the invasion fitness in all possible pairs of mutant and resident architectures (α F,m , α F,r ), with α F,m and α F,r taken among 65 equally spread discrete values in the range [0, 1]. Since the environment state fluctuates cyclically, the equilibrium frequencies of the four combinations of possible (epi)genotype values (X F , X S ) are periodic functions of time (see Section 2.2 and Appendix A). We thus calculated these frequencies at each of 64 timesteps equally spread over a period T . We then calculated, at each timestep and for each pair (α F,m , α F,r ), the growth rate of the mutant. The instantaneous growth rate at time t of the mutant in a given genetic background can be calculated as:

W m,X F ,X S (t) = 1 -s |(1 -α F,m )X S + α F,m X F -Φ * (t)| γ , (2)
with Φ * (t) denoting the optimal phenotype at time t (either Φ * A or Φ * B ). Assuming that the architecture is determined by a modifier locus linked with those controlling the determinants, each combination created at the time of a mutation has its own independent dynamics and long-term relative growth rate, which is not impacted by the change in frequencies in the population and can thus be calculated as the arithmetic mean of the W m,X F ,X S (t) over the period.

The average invasion fitness of a mutant's architecture depends on its probability to appear in each background, represented by background frequencies in the following equation:

W (m, r) =f (0,0) (t)W (m, 0, 0) + f (1,0) (t)W (m, 1, 0) + f (0,1) (t)W (m, 0, 1) + f (1,1) (t)W (m, 1, 1) (3)
We represent the mutants' relative invasion fitnesses for our 65×65 pairs (α m , α r ) in so-called Pairwise Invasibility Plots (PIP), from which we identify singular strategies and characterize their evolutionary stability (according to [START_REF] Geritz | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF].

Simulations in infinite and finite populations

The numerical simulations for infinite populations are based on the discretisation of the continuous model described below. Although the model is valid for any number of inheritance systems, the numerical results are shown only for two inheritance systems in the case of infinite populations. Then, the population is divided in four groups according to their (epi)genotype: g (0,0) , g (1,0) , g (0,1) and g (1,1) . Further, each subpopulation is indexed by individual architecture α F (simply denoted α in the following lines). The density within each subpopulation is driven by

∂ t g X (t, α) = 1 -s|Φ X (α) -Φ * (t)| γ -ρ(t) g X (t, α) + (M G) X (t, α) -µ α g X (t, α) + µ α 1 0 K(α -y, α)g X (t, y)dy. (4)
In this equation:

• 1 -s|Φ X (α) -Φ * (t)| γ -ρ(t)
is the percapita growth rate, including selection and density dependent saturation where ρ(t) is the total population density :

ρ(t) = 1 0 g (0,0) (t, α) + g (1,0) (t, α) + g (0,1) (t, α) + g (1,1) (t, α) dα.
• (M G) X (t, α) describes the (epi)genotype mutation dynamics through the matrix operation M G defined by

M G =          -µ S -µ F µ F µ S 0 µ F -µ S -µ F 0 µ S µ S 0 -µ S -µ F µ F 0 µ S µ F -µ S -µ F                   g (0,0) g (1,0) g (0,1) g (1,1)         
• µ α g X (t, α) represents the mutation toward other trait architecture where µ α is the mutation rate on α,

• µ α 1 0 K(α -y, α)g X (t, y)dy corresponds to the mutation from other trait architecture to the current one; in this term K is the mutation kernel that is assumed to be a Gaussian distribution centered on the progenitors trait architecture, conditioned to take values in [0, 1].

The values of α are uniformly discretised in [0, 1]. The time discretisation is based on an explicit Runge-Kutta (RK2) scheme to reduce computational time and improve accuracy. The integral terms are computed according to standard quadrature methods.

The model ( 4) is non-linear due to the density-dependent saturation. Nevertheless, it can be recast into a linear problem, and, in the periodic setting, the Lyapunov exponent λ(α) coincides with the averaged density over one period 1 T T 0 ρ(t)dt in the case of a monomorphic population when mutations are neglected, see Appendix A.

For finite populations, the stochastic birth and death process described in Section 2.1 is simulated using a Monte-Carlo approach. We opted for a fixed time step discretisation to reduce computational time. Contrary to the infinite population density model (Equation 4), the population remains finite and fluctuates around the carrying capacity K. The generalization of the model to any number of inheritance systems is straightforward, at no computational cost, and so the Monte-Carlo approach is suitable when the number of inheritance systems is larger than three.

Results and discussion

The evolution of the trait architecture α may stem from two distinct effects: (i) its impact on the average fitness across the two environments, regardless of the various mutation rates and (ii) its control of the phenotypic contribution of inheritance systems with different mutation rates and thus presumably different variational properties. Our focus is on the second effect, that is, on the implications of the mutation rates on the trait architecture. We thus begin our analysis (Section 3.1) by setting γ = 1 such that the first effect is removed. Indeed, the averaged fitness does not depend on the phenotype when γ = 1, that is,

in mathematical terms: 1 -s 2 (|Φ -Φ * A | 1 + |Φ -Φ * B | 1 ) = 1 -s 2 (|Φ -0| 1 + |Φ -1| 1 ) = 1 -s 2 for all Φ.
In Section 3.2, we then investigate the sensitivity of our conclusions to changes in the shape of the trait fitness function.

3.1 Architectures evolve in response to the degree of environmental instability

Composite architectures are selected against when the fitness function is linear

To assess which architecture (pure or composite) is optimal under various degrees of environmental instability, we first simulated temporal dynamics in the simple case I = 2 (two inheritance systems, fast and slow, with respective mutation rates µ F and µ S ). In this particular case, we recall that respective

(α) = W (α) max W
is represented as a function of the period T , for three different values of αF = 0, 0.5, 1. The maximal fitness is always attained at extremal values αF = 0 or αF = 1. The analytical value Tc of the transition between 'unstable' and 'stable' environments is reported by a cross-mark (Equation 5). (b) The optimal value (α * F = 0 or α * F = 1) is plotted as a function of the period T for different values of the fast mutation rate µF . Again, the switching time Tc is well approximated by the analytical expression (Equation 5) when µF is sufficiently distant from the values of s (10 -1 ) and µS (10 -4 ).

In fact, the counter-selection of intermediate α F values (and symmetrically, of intermediate α S values) turned out to be a theorem, whose statement is given below:

Assume that the trait fitness function is linear (γ = 1), and I ≥ 2. Then the Lyapunov exponent λ(α) is a convex function of the trait architecture α.

This theorem is a consequence of the convexity properties of the spectral radius with respect to the diagonal coefficients (see e.g. [START_REF] Kingman | A convexity property of positive matrices[END_REF][START_REF] Cohen | Convexity of the dominant eigenvalue of an essentially nonnegative matrix[END_REF]. A proof can be found in Appendix B for the sake of completeness. An immediate consequence of this theorem is that optimal values of the trait architecture are reached at some extremal point of the simplex set {α = (α i ) : (∀i)α i ≥ 0, and i α i = 1}.

Alternatively speaking, composite trait architectures should always be selected against when the selection function is linear.

Another implication of this theorem is the occurrence of sharp transitions, at specific threshold periods, corresponding to switches between optimal trait architectures α * . In the case I = 2, there is one such switch at a critical value T = T c , quantitatively separating unstable environments (T < T c ) where the Lyapunov exponent λ(α) is maximal at α F = 1 (full contribution of the fast inheritance system), from stable ones (T > T c ) where the Lyapunov exponent λ(α) is maximal at α F = 0 (full contribution of the slow inheritance medium), see Figure 3(b). The value of T c can be computed analytically in the regime µ S ≪ µ F , by seeking the value of T such that λ(0) = λ(1). This was performed, and lead to the following expression (after some further simplifications, see Appendix C for details):

T c 2 ≈ 1 µ F log µ F µ S 1 + µ F s (5)
This approximation is valid in the regime µ S ≪ µ F ≪ s, that is, when the mutation rates are far apart, and much below the maximal fitness defect s. It could be further simplified by neglecting µ F s , because selection is assumed to be strong enough so that it has limited influence on the threshold period

T c .
As mentioned above, non-trivial outcomes are expected in the intermediate regime, when the rate of environmental change falls in the range of mutation rates. Formula (5) makes this expectation more precise as the switch occurs when the duration of one stasis (the half-period T 2 ) echoes the fast mutation rate, up to a logarithmic correction. The latter correction is essential as it contains the ratio between the two mutation rates (fast versus slow). Interestingly, although selection is obviously essential for this process, its strength s only weakly influences the value of the switching period (no dependency at the leading order in formula ( 5)). With more than two determinants (differing in their mutation rates), several switching periods are expected, as exemplified in the case I = 3 (see Appendix D, Figure 7).

Epistasis and contingency in the evolution of trait architectures

A mutant genetic architecture may be more or less successful at invading particular (epi)genetic backgrounds. For example, putting more weight on the fast inheritance system may be deleterious if this system is not in the appropriate state in a certain environment (see Section 2). This creates an evolutionary interaction between a resident genetic architecture and the mutant architecture. Importantly, this is ignored in the Lyapunov exponent approach, where the fitnesses of different architectures are compared after each has reached its own (epi)genetic equilibrium frequencies. In contrast, this is well captured by the adaptive dynamics framework. In this framework, the mutant invasion dynamics is approximated through its growth rate at low frequency within a resident population under deterministic dynamics (infinite population size approximation), that is when it is so rare that it does not itself impact the resident's dynamics. This growth rate, calculated as detailed in Section 2, is represented in Figure 4(a) for 4225 (65 2 ) pairs of mutant and resident strategies.

Based on this approach, we identified two types of singular strategies, as represented in Figure 4 q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 period of environmental change, T α q q q q q q q q q q q CSS repellor 10 2 10 3 10 4 (b) starting below or above the repellor will be markedly different. If the initial genotype has an architecture α F below the repellor, e.g. α repellor F ≈ 0.6 in Figure 4(a), mutants with lower α F will tend to invade and reach fixation. Because residents with higher α F consistently fail to invade on average, these evolutionary dynamics will thus lead α F downwards until α F = 0 is reached. The reverse occurs for starting points above the repellor. Thus, the existence of a repellor creates a contingency of the evolutionary dynamics on the initial genotype; whether, or rather how frequently it can be overcome by drift will depend on population features (but see Section 3.3) and the distribution of mutation effects on the trait's architecture.

But it should be noted that selection will consistently act against such crossing of the repellor, as shown by the large white areas in Figure 4(a) for which the probability of fixation is lower than that of a neutral mutant.

We next investigate how the types of singular strategies and their distribution change as a function of environmental instability, see Figure 4(b). When the environment is scarcely changing (up to T ≈ 256), two CSS attractors (at α F = 0 and α F = 1) and one repellor are present and only the position of the repellor changes. In environments that change very rarely, the repellor is very close to α F = 1, such that one may expect the evolution of architectures towards α F = 0, that is, with the lowest possible contribution of the most unstable inheritance system. This conclusion is similar to that reached through the Lyapunov exponent approach, and in line with the results of (Rajon and Charlat 2019). Under higher environmental instability, the repellor moves to intermediate values of α F such that the outcome will depend on the initial genotype. This is again compatible with the results of [START_REF] Rajon | In)exhaustible Suppliers for Evolution? Epistatic Selection Tunes the Adaptive Potential of Nongenetic Inheritance[END_REF] reporting an increase of the average α F in similar contexts, especially since the initial (epi)genotypes were uniformly distributed between 0 and 1 in their study. Decreasing further the period of environmental change below a threshold of T ≈ 256, only the CSS attractor at α F equal or close to 1 remains.

The aforementioned conclusions derived from the adaptive dynamics were supported by Monte-Carlo numerical simulations initialized with monomorphic populations, that is, where all individuals initially have the same trait architecture: either α F = 0.2 as in Figure 5 In contrast with the linear case, non-composite architectures can result in a higher relative fitness ω, especially when the environment is very unstable, or when γ is large. (c) Singular strategies observed under a concave selection function (γ > 1) and different periods of environmental change). Singular strategies are either repellors (empty dots) or CSS attractors of the evolutionary dynamics (plain dots or triangles).

λ(1). Since the intensity of selection does not depend on γ at the extremal architectures (either 1 or 1 -s, and vice-versa), the balance λ(0) = λ(1) does not involve the parameter γ.

It is important to notice that there is no corresponding statement regarding concavity, that is, the trait architecture fitness does not inherit concavity from the trait fitness function in general. For instance, the selection function in the case γ = 1 (Section 3.1.1) is both convex and concave because it is linear.

However, the trait fitness architecture is strictly convex in this case. It is the transition between several (epi)genetic states by mutation, and the environmental variation which generates such strict convexity.

Since the convexity property is a strong feature, which holds true even in the case of a linear selection function, we expect that convexity of the Lyapunov exponent persists beyond the linear regime, even if the trait fitness function is not convex anymore, say slightly concave. Robustness of the previous conclusions with respect to the shape of the fitness function is the purpose of the next section.

Concave fitness functions may select for composite architectures

Contrary to convexity, concave fitness functions may result in positive selection for a balanced contribution of various determinants -that is 0 < α F < 1 -compare Figure 6 formula for the critical value T c in the concave case, we see clearly that it does not depend on the shape parameter (γ). Although this may appear surprising at first glance, it might be expected from the fact that the critical value does not depend on the shape parameter (γ), nor on the intensity of selection (s)

at the leading order in the convex case, see Equation ( 5).

Notably, non-composite architectures can still be selected for, for moderate shape parameter (γ ≲ 1.5) and period T close enough to the switching period T c . For large γ, the fastest inheritance system (α * F = 1) is never selected for, and the slowest one is selected for only in very stable environments (T ≫ T c ).

Regarding the adaptive dynamics, the results in Figure 6

Sensitivity to population size: a new case of sign inversion

We used simulations to assess the predictions of the adaptive dynamics approach under more or less effective selective regimes, that is, under various populations sizes. Simulations behave as expected when the population size is sufficiently large (K = 4096), that is α F goes to 0 when the population is initiated with α F below the repellor and to 1 otherwise. In very small populations, selection is overall inefficient and α F tends to become uniformly distributed, independently of its initial value (K = 32 in Figure 5). This is also expected since genetic drift tends to dominate over selection in these conditions. But in populations of intermediate sizes (K = 512 in Figure 5), simulations reveal a complex pattern. When the initial value α F = 0.8 is above the repellor, it evolves toward 1, as expected. Yet, surprisingly, α F also evolves toward 1 when it is initiated below the repellor, at initial value α F = 0.2. In other words, we found that not only the efficiency of selection, but also its direction, depends on the population size (compare panels (d) and (e) of Figure 5, where the only parameter that differs is K).

We reached the following interpretation for this surprising pattern. With this explanation in mind, we realized that other authors have recently shed light on a very similar process that they coined 'sign inversion' [START_REF] Raynes | Sign of selection on mutation rate modifiers depends on population size[END_REF][START_REF] Raynes | When good mutations go bad: how population size can change the direction of natural selection[END_REF].

This expression appropriately emphasizes the main point of the explanation: changes in population size can not only make selection more or less efficient, but may also inverse its direction in systems where two opposing pressures are at play and differ in their intensity and duration. In large populations, both pressures are effective, but weak and long lasting selection may dominate; under a certain threshold population size, only the strong but infrequent selection is effective. Following the discovery of this process in a study focusing on the evolution of mutation rates, [START_REF] Raynes | Sign of selection on mutation rate modifiers depends on population size[END_REF] have noted that this process may be far more than anecdotal, because many traits can be subject to opposing selective pressures, for example because of temporal variations in the environment. Our study provides one more example of situations where this process would take place.

Conclusion

Our analysis was aimed at assessing the hypothesis that heterogeneity in the mutation rates of a trait's determinants (of whatever nature) may affect the evolution of its architecture, that is, of the respective contributions of the different determinants. Our main conclusion is that in a variety of conditions, in both stable and unstable environments, non-composite architectures are effectively selected for, meaning that trait values would tend to depend upon a single determinant, or several determinants characterized by similar mutation rates. This sheds light on prior simulation results [START_REF] Rajon | In)exhaustible Suppliers for Evolution? Epistatic Selection Tunes the Adaptive Potential of Nongenetic Inheritance[END_REF] as we now understand how the degree of environmental variation conditions which determinant, rather than which mixture of determinants, takes part to the trait architecture. This focus on the mutational portfolio of traits contributes to the general understanding of the evolution of a trait's architecture, as a complement to the more common focus on the number of contributing genes, their epistatic and pleiotropic properties, and the distribution of their respective effects.

Let us emphasize that a non-composite trait architecture may still be "complex" in the sense of involving many determinants (e.g. [START_REF] Ungerer | Quantitative trait loci for inflorescence development in Arabidopsis thaliana[END_REF][START_REF] Flint | Genetic architecture of quantitative traits in mice, flies, and humans[END_REF][START_REF] Kemper | Genetic architecture of body size in mammals[END_REF][START_REF] Yengo | A saturated map of common genetic variants associated with human height[END_REF]. For example, Rajon and Plotkin 2013, in a model assuming a single mutation rate, have shown that many genes should contribute to traits subject to a stabilizing selection of mild intensity. In such a context, where the environment is kept stable, our model predicts that the architecture should be non-composite. Overall, a complex yet non composite architecture would then be expected, involving the most stable determinants.

The impact of environmental heterogeneity on the evolution of trait architecture has also been the subject of previous studies. Yeaman and Whitlock 2011 have shown that spatial differences in traits optima leads the evolution of recombination landscapes toward simplified architecture where clusters of genes have large impact of the phenotype. Using yet another framework, the results of [START_REF] Hansen | Evolution of Genetic Architecture under Directional Selection[END_REF]) suggest that directional selection, producing an ever going requirement for adaptation should result in an increase in mutational effects through changes in epistatic interactions, provided that mutations are rare enough.

In our study, we focused on temporal environmental fluctuations on the evolution of the trait architecture. More specifically, we assumed that environmental changes reverse the direction of selection, with fitness gains modeled as concave, linear, or strictly convex functions. These shapes were found critical: only a concave fitness function resulted in composite architectures, provided that the degree of instability overcomes a threshold. The same conclusion holds for a mixture of concave and convex fitness shapes, like the typical case of a Gaussian function centered on a moving optimum, see Appendix D, Figure 8(a).

In such a case, the amplitude of the fluctuations may be critical because Gaussian functions are concave close to the center, but convex at the border, as can be seen in Appendix D, Figure 8(b).

While the shape of fitness functions acting on the trait is a key parameter, it will only be relevant if selection on the architecture is effective, which ultimately depends on its strength, in relation with population size. We developed complementary approaches to investigate the strength of selection on the trait architecture. The Lyapunov exponent approach revealed conditions where composite architecture are not optimal, but it missed the notion that the evolution of trait architecture is driven by the background variations, accumulated by different determinants through their distinct mutation rates and across environmental fluctuations. The adaptive dynamics approach captured this effect, by assessing the average strength of selection over an environmental fluctuation period. But Monte Carlo simulations in populations of finite size revealed a singular aspect of these dynamics, related to the effect of drift. Simulations in large populations were in agreement with the adaptive dynamics predictions, and those in very small populations were as expected driven by drift only. Yet, in populations of intermediate sizes, we found that the trait architecture evolves deterministically in the opposite direction as in large populations. We interpret this result as a new occurrence of 'sign inversion' [START_REF] Raynes | Sign of selection on mutation rate modifiers depends on population size[END_REF][START_REF] Raynes | When good mutations go bad: how population size can change the direction of natural selection[END_REF] Weinreich 2021), stemming from the periodic alternation between long time windows of weak selection for a heavy contribution of stable inheritance systems (inefficient in populations of intermediate sizes) and short time windows of strong selection for a heavy contribution of unstable systems. This result further complexifies our understanding of the evolution of genetic architectures, arguing in favor of a systematic use of population genetics models and, as an aside, for further mathematical analyses that would clarify and generalize the conditions for the occurrence of sign inversion.

Appendices A Model reduction for Lyapunov exponent computation

In this section, we aim to detail the links between the continuous model presented in Section 2.4 and formalized by Equation ( 4) and the Lyapunov exponent computation. Indeed, the dynamic is driven by the leading eigenvalues in a linear ordinary differential equation system. The presence of the integral term for the trait architecture mutation and the logistic growth rate to model the density dependent saturation make our model non-linear. However, in the case of a monomorphic population, when the trait architecture mutations are neglected, we demonstrate a correspondance with the Floquet spectral theory of periodic linear systems, such that the average of the density over one period, 1 T T 0 ρ(t)dt, and the Lyapunov exponent defined in Equation ( 1), coincide.

Assume that mutations on the trait architecture α are negligible. For a given trait architecture α and group X ∈ (0, 0), (1, 0), (0, 1), (1, 1) the equation ( 4) becomes

∂ t g X (t, α) = 1 -s|Φ X (α) -Φ * (t)| γ -ρ(t) g X (t, α) + (M G) X (t, α). (6)
Depending on the optimal phenotype Φ * , the system gathering equations (6) writes

       ∂ t G(t, α) = (A -ρ(t)I)G(t, α), if Φ * (t) = 0, ∂ t G(t, α) = (B -ρ(t)I)G(t, α), if Φ * (t) = 1. (7)
Assume that both matrix A and B can be diagonalised. Thus, there exist a constant invertible matrix P A (resp. P B ) and a diagonal matrix D A (resp. D B ) such that

A = P A D A P -1 A (resp. B = P B D B P -1 B ).
So the systems (7) becomes

       ∂ t P -1 A G(t, α) = D A -ρ(t)I P -1 A G(t, α), if Φ * (t) = 0, ∂ t P -1 B G(t, α) = D B -ρ(t)I P -1 B G(t, α), if Φ * (t) = 1. ( 8
)
With the notations H A (t, α) = P -1 A G(t, α) and H B (t, α) = P -1 B G(t, α) this systems writes

       ∂ t (H A (t, α)) = D A -ρ(t)I H A (t, α), if Φ * (t) = 0, ∂ t (H B (t, α)) = D B -ρ(t)I H B (t, α), if Φ * (t) = 1. (9)
Relations between H components: Let H A,i denotes the i th components of H A and d A,i the i th diagonal term of D A . Assume that the H A,i never vanishes, then

∂ t H A,i (t, α) H A,i (t, α) -d A,i = ρ(t). (10)
Thus, for any pair (i, j)

∂ t H A,i (t, α) H A,i (t, α) -d A,i = ∂ t H A,j (t, α) H A,j (t, α) -d A,j (11)
=⇒ H A,i (t, α) = H A,j (t, α) H 0 A,i H 0 A,j e (dA,i-dA,j)(t-t0) . ( 12
)
with H 0 A = H A (t 0 , α). Therefore, the relative dynamics of the components of the vector H depend on the difference between the eigenvalues of A. By the same reasoning and using similar notations it also comes

H B,i (t, α) = H B,j (t, α) H 0 B,i H 0 B,j
e (dB,i-dB,j)(t-t0) .

(13) Solution characterization using total population size: According to systems (9), as long as Φ * (t) remains constant on the time interval [t 0 , t] with t > t 0 the solutions are

       H A (t, α) = exp D A (t -t 0 ) - t t 0 ρ(s)ds H 0 A , if Φ * (t) = 0, H B (t, α) = exp D B (t -t 0 ) - t t 0 ρ(s)ds H 0 B , if Φ * (t) = 1 (14)
with H 0 A = H A (t 0 , α) and H 0 B = H B (t 0 , α). Note that these formulas rely on

       G(t, α) = P A exp D A (t -t 0 ) - t t 0 ρ(s)ds P -1 A G 0 , if Φ * (t) = 0, G(t, α) = P B exp D B (t -t 0 ) - t t 0 ρ(s)ds P -1 B G 0 , if Φ * (t) = 1. ( 15
)
where G 0 = (t 0 , α). These formulas simply into

       G(t, α) = exp - t t 0 ρ(s)ds e A(t-t 0 ) G 0 , if Φ * (t) = 0, G(t, α) = exp - t t 0 ρ(s)ds e B(t-t 0 ) G 0 , if Φ * (t) = 1. (16)
A.1 Groups dynamic in stable environment

Without loss of generality, let's consider the case of Φ * (t) = 0 and drop the subscript A to shorten notations. Using G(t, α) = P H(t, α) one can deduce that the expression of G i the i th components of G is

G i (t, α) = n j=1 P i,j H j (t, α) = H ℓ (t, α) H 0 ℓ n j=1 P i,j H 0 j e (d j -d ℓ )(t-t 0 ) . ( 17
)
where ℓ is an arbitrary index to be fixed later. Thus, the total population size is

ρ(t) = n i=1 G i (t, α) = H ℓ (t, α) H 0 ℓ n i=1 n j=1 P i,j H 0 j e (d j -d ℓ )(t-t 0 ) . (18)
Therefore, the relative group size within the total population size is

G i (t, α) ρ(t) = n j=1 P i,j H 0 j e (d j -d ℓ )(t-t 0 ) n i=1 n j=1 P i,j H 0 j e (d j -d ℓ )(t-t 0 ) (19)
= n j=1 P i,j H 0 j e d j (t-t 0 ) n i=1 n j=1 P i,j H 0 j e d j (t-t 0 )

(20)

Thus, the relative group size dynamic is fully determined by the eigenvalues of matrix A. Moreover, this shows that the whole dynamic is driven by the A matrix leading eigenvalue:

G i (t, α) ρ(t) ≃ P i,ℓ H 0 ℓ e d ℓ (t-t 0 ) n i=1 P i,ℓ H 0 ℓ e d ℓ (t-t 0 ) = P i,ℓ n i=1 P i,ℓ (21)
assuming that d ℓ is the leading eigenvalue. Note that P i,ℓ n i=1 P i,ℓ correspond to the normalised value of the i th component of the eigenvector associated to d ℓ . That is, the normalised eigenvector associated with the largest eigenvalue of A approximates the relative group sizes.

A.2 Groups dynamic in periodic environment

Let's consider a periodic environment of period T . Without loose of generality let assume that Φ * (t) = 0 for t ∈ 0, T 2 and Φ * (t) = 1 for t ∈ T 2 , T . Thus, starting from G 0 = G(0, α) and using formulas ( 16) the Up to renormalization, the eigenvector G 0 can be interpreted as the relative group sizes within the whole population. Indeed, the relative size of a group ℓ is simply

solution at t = T is G(T, α) = exp -
(G 0 ) ℓ i (G 0 ) i .

B Proof of convexity

This section is devoted to the mathematical proof of our main convexity result, that is the convexity of the Lyapunov exponent λ(α) in the case of a convex selection function (γ ≤ 1). This result can be recast in the framework of convexity of dominant eigenvalues of non-negativity matrices. To apply such theory, one should be careful about the convex combination of coefficients during convex interpolation.

Indeed, the diagonal coefficients and the off-diagonal coefficients do not play the same role, for as they are not subject to the same constraints (non-negativity constraints off the diagonal). More precisely, convex combinations are arithmetical on the diagonal, but geometrical off the diagonal (see the θ-interpolation below). The situation is facilitated in our model as the selection acts only on the diagonal terms via a trait-dependent mortality rate. Nevertheless, we shall present a more general result to emphasize this discrepancy between diagonal and off-diagonal coefficients.

We claim that the following Theorem encompasses our main convexity result:

Let A, B two matrices with positive entries off the diagonal1 , that is

∀(i, j) a ij > 0, b ij > 0.
Let θ ∈ (0, 1). Consider the matrix C θ defined by the following term-by-term interpolation:

       c ii = (1 -θ)a ii + θb ii c ij = a 1-θ ij b θ ij (i ̸ = j)
Then the dominant (Perron) eigenvalues λ(A), λ(B), λ(C θ ) satisfy the following convex inequality:

λ(C θ ) ≤ (1 -θ)λ(A) + θλ(B) (22)
We make the following comments:

• This result is classically credited to [START_REF] Kingman | A convexity property of positive matrices[END_REF]).

• The result can be extended straightforwardly to irreducible matrices by a limiting argument.

• We are interested in the particular case where A and B differ only by their diagonal entries, that is, a ij = b ij if i ̸ = j. In this case, the interpolation c ij = a 1-θ ij b θ ij is trivial, and the convex inequality ( 22) can be reformulated as: λ((1 -θ)A + θB) ≤ (1 -θ)λ(A) + θλ(B). Alternatively speaking, the dominant eigenvalue is convex with respect to its diagonal.

• In fact, the proof can be extended to the Floquet eigenvalues in the periodic setting, which is precisely our focus:

Let A(t), B(t) a pair of time-dependent, periodic, matrices with positive entries off the diagonal. Let θ ∈ (0, 1). Consider the matrix C θ (t) defined pointwise as follows:

       c ii (t) = (1 -θ)a ii (t) + θb ii (t) c ij (t) = a 1-θ ij (t)b θ ij (t) (i ̸ = j)
Then the dominant (Floquet) eigenvalues λ(A), λ(B), λ(C θ ) satisfy the following convex inequality:

λ(C θ ) ≤ (1 -θ)λ(A) + θλ(B) (23)
We present first the proof of the static case (Perron eigenvalues), then we extend the method to the time-periodic setting (Floquet eigenvalues).

Our approach is based on [START_REF] Clairambault | Circadian rhythm and cell population growth[END_REF], see also [START_REF] Cohen | Convexity of the dominant eigenvalue of an essentially nonnegative matrix[END_REF]).

The static case

We need the following eigen-elements: X (resp. Y ) the positive (right) eigenvector of A (resp. B), and Φ the positive (left) eigenvector of C. We have accordingly:

       AX = λ(A)X BY = λ(B)Y , Φ T C = λ(C)Φ T (24)
We define Z by the pointwise interpolation:

z i = x 1-θ i y θ i .
Then, we evaluate the Collatz-Wielandt ratio:

j c ij z j z i = (1 -θ)a ii + θb ii + j̸ =i a 1-θ ij b θ ij x 1-θ j y θ j x 1-θ i y θ i ≤ (1 -θ)a ii + θb ii + (1 -θ)   j̸ =i a ij x j x i   + θ   j̸ =i b ij y j y i   = (1 -θ)λ(A) + θλ(B)
where we have used the Hölder inequality in the second line, and the fact that X and Y are eigenvectors in the last line. We deduce that

CZ ≤ ((1 -θ)λ(A) + θλ(B))Z
Finally, multiplying by Φ T , we find:

λ(C)Φ T Z ≤ ((1 -θ)λ(A) + θλ(B))Φ T Z, (25) hence
, λ(C) ≤ (1 -θ)λ(A) + θλ(B).

The time-periodic case

The proof is almost identical, but the fact that we have to consider time-periodic eigenvectors X(t) (resp.

Y (t)) and Φ(t) such that:

       -Ẋ(t) + A(t)X(t) = λ(A)X(t) -Ẏ (t) + B(t)Y (t) = λ(B)Y (t) , ΦT (t) + Φ T (t)C(t) = λ(C)Φ T (t) (26)
We define again Z(t) pointwise: z i (t) = x 1-θ i (t)y θ i (t). Then, we evaluate the Collatz-Wielandt ratio:

- żi z i + j c ij z j z i = -(1 -θ) ẋi x i -θ ẏi y i + (1 -θ)a ii + θb ii + j̸ =i a 1-θ ij b θ ij x 1-θ j y θ j x 1-θ i y θ i ≤ (1 -θ)   λ(A) - j a ij x j x i   + θ   λ(B) - j b ij y j y i   + (1 -θ)a ii + θb ii + (1 -θ)   j̸ =i a ij x j x i   + θ   j̸ =i b ij y j y i   = (1 -θ)λ(A) + θλ(B)
We deduce that -Ż(t) + C(t)Z(t) ≤ ((1 -θ)λ(A) + θλ(B))Z(t)

To conclude, we multiply by ϕ T and integrate other one period:

- T 0 -Φ T (t) Ż(t) dt + T 0 Φ T (t)C(t)Z(t) dt ≤ ((1 -θ)λ(A) + θλ(B)) T 0 Φ T (t)Z(t) dt
Integrating by parts the first term (using periodicity), we find

C Computation of the critical time T c

This section presents the analytical computation of the critical time T c of the transition between 'unstable' and 'stable'. According to our main results, see Section 3.1.1: when the trait fitness function is linear (γ = 1), then the Lyapunov exponent λ(α) is a convex function of the trait architecture α. Consequently, for the two determinant case the optimal trait architecture is either2 α F = 0, or α F = 1 and the transition occurs when the Lyapunov exponents satisfy λ(0) = λ(1). Thus, the critical time T c can be computed by solving this equation. However, according to Equation (1), the Lyapunov exponent is the largest eigenvalue of a product of matrix exponential. Thus, its computation is very tedious so we relied on Maple (2022) to perform the preliminary steps of the calculation.

As mentioned above, the goal is to solve λ(0) = λ(1) thus the first step consists in determining the Lyapunov exponent for α F = 0 and α F = 1. To shorten the notation, we abbreviate α F to α. According to Maple, the eigenvalues of exp The mutation rate µ F is strictly positive so e -4µ F < 1 and thus the largest eigenvalue is either Λ(0) 1 or Λ(0) 2 . Besides, ξ S simplifies into ξ S = µ S e -τ S 2 -1 e -τ S 2 + 1 2 µ 2 S + s 2 e -τ S 2 > 0.

Thus, when α = 0 the largest eigenvalue is λ(0) T = Λ(0) 1 = exp T 2 (2 -2µ S -s) + τ S 2 4µ 2 S + s 2 2µ 2 S e -τ S + s 2 e -τ S 2 + 2µ 2 S + 2ξ S .

Similarly, using Maple, we found that the eigenvalues of exp AT 2 exp BT 2 when α = 1 are

Λ(1) = exp T 2 (2 -2µ F -s) + τ F 2 4µ 2 F + s 2          2µ 2 F e -τ F + s 2 e -τ F 2 + 2µ 2 F + 2ξ F 2µ 2 F e -τ F + s 2 e -τ F 2 + 2µ 2 F -2ξ F 2µ 2 F e -τ F + s 2 e -τ F 2 + 2µ 2 F + 2ξ F e -4µ F 2µ 2 F e -τ F + s 2 e -τ F 2 + 2µ 2 F -2ξ F e -4µ F          with τ F = T 4µ 2 F + s 2 ,
ξ F = µ 2 F e -2τ S µ 2 F -2 e -τ S µ 2 F -2 e -τ S s 2 + e -τ S 2 s 2 + e -3τ S 2 s 2 + µ 2 F .

As previously ξ F simplifies into ξ F = µ F e -τ S 2 -1 e -τ S 2 + 1 2 µ 2 F + s 2 e -τ S 2 > 0 and using the same arguments we found that when α = 1 the largest eigenvalue is

λ(1) T = Λ(1) 1 = exp T 2 (2 -2µ F -s) + τ F 2 4µ 2 F + s 2 2µ 2 F e -τ F + s 2 e -τ F 2 + 2µ 2 F + 2ξ F .
Therefore, the critical time T c can be computed by solving Λ(0) 1 = Λ(1) 1 . In practice, we were unable to solve this equation. Even the numerical approximation of the solution requires particular care due to the presence of large terms in the exponentials. Nevertheless, the terms of the form e -cτ i with c > 0 a constant and i ∈ {F, S} are negligible in the regime µ F ≪ s. Thus, for α = 0 it comes ξ S ≈ µ 2 S and then

Λ(0) 1 ≈ 4µ 2 S 4µ 2 S + s 2 exp T 2 2 -2µ S -s + τ S 2 .
Similarly, for α = 0 it comes ξ S ≈ µ 2 F and then

Λ(1) 1 ≈ 4µ 2 F 4µ 2 F + s 2 exp T 2 2 -2µ F -s + τ F 2 .
So, the approximated critical time T c is the solution of the equation:

4µ 2 S 4µ 2 S + s 2 exp T c 2 2 -2µ S -s + τ S 2 = 4µ 2 F 4µ 2 F + s 2 exp T c 2 2 -2µ F -s + τ F 2 ⇐⇒ exp T c 2 2 -2µ S -s + τ S 2 - T c 2 2 -2µ F -s - τ F 2 = 4µ 2 F 4µ 2 S + s 2 4µ 2 S 4µ 2 F + s 2 ⇐⇒ exp 1 2 (τ S -τ F ) + T c (µ F -µ S ) = 4µ 2 F 4µ 2 S + s 2 4µ 2 S 4µ 2 F + s 2 ⇐⇒ 1 2 T c 4µ 2 S + s 2 -T c 4µ 2 F + s 2 + T c (µ F -µ S ) = log 4µ 2 F 4µ 2 S + s 2 4µ 2 S 4µ 2 F + s 2 .
Thus, the critical time is approximately

T c ≈ 2 log 4µ 2 F (4µ 2 S +s 2 ) 4µ 2 S (4µ 2 F +s 2 ) 4µ 2 S + s 2 -4µ 2 F + s 2 + 2(µ F -µ S )
.

However, as it stands, the contribution of each parameter to the critical time remains complex. Thus, according to Taylor expansions in the regime µ S ≪ µ F ≪ s, this formula simplifies into

T c 2 ≈ 1 µ F log µ F µ S 1 + µ F s .

Figure 1 :

 1 Figure1: Shape of the fitness function 1 -s|Φ -Φ * | γ , for γ = 1 (plain green, linear case), γ = 1/2 (dashed plain red, strictly convex case) and γ = 2 (dotted blue, strictly concave case).

Figure 2 :Figure 3 :

 23 Figure 2: Temporal dynamics of the trait architecture with polymorphic initial configuration. (a) Simulation of the deterministic PDE model (Equation4), corresponding to an infinite population size, with a polymorphic initial configuration where all αF values have the same frequency. The mean αF in the population is plotted over time for two values of the time period of environmental change T . We observed the selection towards αF = 1 for the smaller period T = 900 and the selection towards αF = 0 for the larger period T = 1200. (b-c) Monte-Carlo simulations associated with finite population sizes, respectively for T = 900 (b) and T = 1200 (c). The equilibrium population size is set to K = 2 16 ≈ 65000. The initial value of αF is drawn uniformly within the population. Each black line represents the trajectory of the averaged value αF within the population over time. The colored lines represent the mean and median of 64 simulation replicates.

  (a): evolutionary repellors from which evolutionary trajectories tend to move away, and convergent stable strategies (CSS) that represent possible long-term evolutionary attractors. Evolutionary trajectories

Figure 4 :

 4 Figure4: Adaptive dynamics of trait architectures. (a) A pairwise invasibility plot showing how evolutionary trajectories of the genetic architecture may depend on historical contingency. Black squares represent mutant/resident pairs where the relative growth rate of the mutant λ(αF m, αF r ) is close to 1 (±10 -9 ), while grey and white squares represent pairs where λ(αF m, αF r ) > 1 and < 1, respectively. A first example evolutionary trajectory is shown in orange: a mutant with a lower contribution of the 'fast' component αF m appears in a population with αF r close to 0.6. The mutant increases in frequency (i.e. it invades; grey square). It is expected to reach fixation (i.e. replace the resident) because the latter cannot invade back when it is rare (white square at the opposite coordinate, red cross). The evolutionary trajectory continues downwards until the architecture αF = 0 is fixed. Starting from a slightly different architecture (αF r ≈ 0.7, in green), the evolutionary trajectory goes in the opposite direction towards αF = 1, showing how the evolution of the genetic architecture may be contingent on the ancestral architecture of a population. Mutants increasing αF in the orange trajectory or decreasing αF in the blue one, not represented, may have appeared and failed to invade. The evolutionary repellor at αF ≈ 0.65 is represented by an empty dot and the two CSS attractors at αF = 0 and 1 by plain dots. Parameter values: γ = 1, T = 2048, s = 0.1. (b) Singular strategies observed under a linear trade-off (γ = 1) and different frequencies of environmental change (x-axis in log scale). Singular strategies are either repellors (empty dots) or CSS attractors of the evolutionary dynamics (plain dots).

Figure 5 :

 5 Figure1 (γ < 1). Indeed, the trait architecture fitness inherits convexity from the trait fitness function, as stated in the following theorem.

Figure 6 :

 6 Figure 6: Main outcomes for concave fitness functions. (a) Temporal dynamics for the deterministic PDE models of infinite population size for three values of the time period T = 100, T = 900, and T = 1200, and γ = 3 (analogous to Figure 2(a)). The mean value of αF in the population is plotted over time. In contrast with the linear case, we observed the selection for non-composite architectures, especially when the environment is very unstable (T = 100). (b) The optimal value α * F is represented as a function of the shape parameter γ > 1 (analogous to Figure 3(b)).In contrast with the linear case, non-composite architectures can result in a higher relative fitness ω, especially when the environment is very unstable, or when γ is large. (c) Singular strategies observed under a concave selection function (γ > 1) and different periods of environmental change). Singular strategies are either repellors (empty dots) or CSS attractors of the evolutionary dynamics (plain dots or triangles).

  (a) with Figure 2(a). Thus, in order to assess the robustness of our main conclusions from Section 3.1.1, we investigated the optimal trait architecture maximizing the Lyapunov exponent λ(α) in the case of a concave selection function (γ ≥ 1), see Figure 6(b). Interestingly, the concept of a critical switching value T c still exists when composite trait architectures are selected for, even for large values of γ. Fast inheritance systems have more weight in the selected architecture (α * F > 0.5) when T < T c . In contrast, slow inheritance systems have more weight in the selected architecture (α * F < 0.5), when T > T c . Moreover, the transition is sharp at T = T c (switch between a clear balance on the side of fast inheritance towards a clear balance on the opposite side), a feature shared with the convex case, compare Figure 6(b) with Figure 3(b). Although we lack a

  (c) are consistent with the previous results.When the selection function is moderately concave (here, γ = 1.5), and the environment changes infrequently (T > 8192), a single CSS attractor is present at α F close to 0. For intermediate regimes of environmental variation, a repellor and the two attractors at respectively low and high α F are present, similar to the case where γ = 1 (compare with Figure6(b)). As environmental variations become more frequent, at T < 1024 a single CSS remains and the expected selected architecture tends to intermediate values of α F , ending up at α F = 0.5 for very frequent changes. As concavity increases to γ = 3, the transition between an architecture with a large contribution of the determinant with the lowest mutation rate (low α F ) in scarcely varying environments to an architecture with a higher α F no longer includes a repellor, nor a very large contribution of the least stable determinant -i.e. α F does not reach values close to 1.

  Assume that the solution becomes periodic, namely G(T, α) = G 0 then the previous formula lead to e the leading eigenvalue is real.

T 0 ( 0 Φ

 00 ΦT (t) + Φ T (t)C(t))Z(t) dt ≤ ((1 -θ)λ(A) + θλ(B)) )Φ T (t)Z(t) dt ≤ ((1 -θ)λ(A) + θλ(B)) T T (t)Z(t) dt Hence, λ(C) ≤ (1 -θ)λ(A) + θλ(B).

  The evolution of α F in fluctuating environments is dynamically driven by two opposite selective pressures: (1) upon shifts to a new environment, high values of α F are selected for, because these give more weight to the "fast" inheritance system, that has a higher chance of being in a state that matches the new environment, or to reach this state through an (epi)mutation; (2) during periods of environmental stasis, low values of α F are selected for, because these reduce the cost of frequent deleterious (epi)mutations. Depending on the frequency of environmental change, one of these processes tends to dominate, provided that populations are large enough. A plausible interpretation for the observed pattern (that is, for a deterministic evolution of α

F toward high values in populations of "intermediate size") relies on noting that the strength and duration of the selection associated with the two processes differ. As a result, one process, subject to stronger selection than the other, may become effective in populations where the other is not. More specifically, what we see here is in our view stemming for a strong (but infrequent) selection toward high α F values upon environmental change and a weak (but durable) selection toward low α F values during periods of environmental stasis. Both are effective in large populations, but weak and durable selection toward low α F dominates. In contrast, in populations of intermediate size, only the strong albeit infrequent selection toward high α F operates.

  Using the relation G(t, α) = PA H A (t, α) and G(t, α) = P B H B (t, α) it comes

	t t 0 ρ(s)ds which is not
	explicitly known.

  S e -τ S + s 2 e -τ S 2 + 2µ 2 S + 2ξ S 2µ 2 S e -τ S + s 2 e -τ S 2 + 2µ 2 S -2ξ S 2µ 2 S e -τ S + s 2 e -τ S 2 + 2µ 2 S + 2ξ S e -4µ F 2µ 2 S e -τ S + s 2 e -τ S 2 + 2µ 2 S -2ξ S e -4µ F S = µ 2 S e -2τ S µ 2 S -2 e -τ S µ 2 S -2 e -τ S s 2 + e -τ S 2 s 2 + e -3τ S 2 s 2 + µ 2 S .

		AT 2	exp BT 2	when α = 0 are
					2µ 2	
	Λ(0) =	exp T 2 (2 -2µ S -s) + τ S 2 4µ 2 S + s 2	       	       
	with
	τ S = T 4µ 2 S + s 2 ,
	ξ

usually referred to as a Metzler matrix

We recall that in the two determinant case α = (αF , αS) with αS = 1 -αF .

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 639638 and grant agreement No 865711). Computations have been performed on the supercomputer facilities of the Mésocentre de calcul de Franche-Comté.

D Supplementary figures

Relative fitness (0/0/1) (0/0.5/0.5) (0/1/0) (0.5/0/0.5) (0.5/0.5/0) (1/0/0) Predicted Tc