Somatic maintenance impacts the evolution of mutation rate

Background The evolution of multi-cellular animals has produced a conspicuous trend toward increased body size. This trend has introduced at least two novel problems: an expected elevated risk of somatic disorders, such as cancer, and declining evolvability due to generally reduced population size, lower reproduction rate and extended generation time. Low population size is widely recognized to explain the high mutation rates in animals by limiting the presumed universally negative selection acting on mutation rates. Results Here, we present evidence from stochastic modeling that the direction and strength of selection acting on mutation rates is highly dependent on the evolution of somatic maintenance, and thus longevity, which modulates the cost of somatic mutations. Conclusions We argue that the impact of the evolution of longevity on mutation rates may have been critical in facilitating animal evolution. Electronic supplementary material The online version of this article (10.1186/s12862-019-1496-y) contains supplementary material, which is available to authorized users.


Background
Increasing body size has been one of the major trends in animal evolution across many taxa, as formulated in Cope's rule [1,2]. The evolution of larger bodies, particularly for land animals and especially for mammals, introduces some fundamentally new evolutionary challenges. The carrying capacity of ecosystems limits biomass per group/species, so larger body size leads to reduced population size. Furthermore, large animals generally demonstrate lower numbers of progeny per time period, longer generation times, and more advanced ages of first reproduction [3,4]. In aggregate, such changes weaken selection that can act on a population and thus should negatively affect evolvability, defined as "an organism's capacity to generate heritable phenotypic variation" [5]. This general reduction in evolvability should, however, be at least partially alleviated by diversity facilitated by sexual reproduction.
The mutation rate (MR) is another critical evolvability parameter [6]. It is believed that selection generally acts to lower MR [7][8][9], and the significantly higher MRs observed in animals compared to unicellular organisms have been argued to result from the reduced power of selection imposed by small population sizes [10][11][12]. As shown for human males, the germline mutation rate has also been proposed to be related to age at reproduction, given continued replication of germ cells [13,14]. An increasing body of evidence indicates that germline (gMR) and somatic (sMR) mutation rates are linked, as they employ the same basic DNA replication and repair machinery, and loss of function mutations in factors important for various DNA repair pathways result in increases in mutation rates for both somatic and germline cells [15][16][17][18][19][20][21] (here we will understand gMR as genetic variation passed over generations, since the mutation rate in germline cells per cell division should be understood as a somatic mutation rate). In particular, Pothof et al. identified 61 genes whose mutation affects sMR in C. elegans, and mutation of 23 of 24 of the genes tested showed similar impacts on gMR [15]. Moreover, Chen et al. examined the effects of transcriptional level of a gene, GC content, histone marks, and replication timing on mutation rate, and while differences were evident (related to transcriptional level), these parameters largely impacted mutation rates similarly in the germline and soma [18]. It should be noted, however, that linear correlations between gMR and sMR are not expected, as some factors, such as the above mentioned differences related to how transcriptional level of a gene impacts transcription-coupled repair [18], should exhibit different effects on gMR and sMR. Some studies also report differences in the absolute values of the mutation rates in the germline and the soma [22], noting that estimates of somatic mutation rates in vivo are inaccurate due to the need to roughly estimate the numbers of underlying cell divisions [23]. While elevated gMR improves evolvability (responsiveness to selection), the ensuing higher sMR should elevate the risk of somatic disorders, such as cancer [24]. For cancer, increasing body size is expected to increase the frequency of oncogenic mutations by increasing the number of target cells [25]. Somatic mutations also contribute to aging and a variety of aging-related diseases [26]. The increased cost of sMR should thus exert negative selective pressure on gMR in larger animals [27,28].
Recent evidence demonstrates that the sMR in some animal tissues can be higher than the rate inferred from observed mutations, because somatic purifying selection is reasonably effective at eliminating damaged somatic cells with substantial genomic damage [29]. Many mechanisms, such as various tumor suppressor gene functions (including DNA damage induced apoptosis) [30], autophagy [31], purifying somatic selection [29,32], and immune surveillance [33], should buffer the costs of somatic mutations and in aggregate promote lifespan extension by maintaining tissue integrity. We will collectively call these mechanismsthe somatic maintenance program (SMP).
We present theoretical evidence from Monte Carlo modeling indicating that somatic maintenance not only improves individuals' survival for large animals by reducing sMR costs, but should have played a crucial role in animal evolution by impacting the evolution of gMR. We show that positive selection for increased body size promotes positive selection for extended longevity by improving SMP. Our results also indicate that positive selection acting on traits that do not impact somatic risks also promotes selection for an improved SMP. In both cases, evolution toward increased gMR was observed because of the reduced sMR cost, which dramatically improved evolvability of the simulated population.

Theoretical introduction to the modeling
We built a stochastic model of evolution in animal populations, incorporating reproduction and survival ( Fig. 1), whereby each individual's trait is inherited with variance proportional to gMR (for code, see Additional file 1: Section 1a). Traits are assumed to be polygenic and exhibit phenotypic variation in the population. In particular, MR is assumed to be a highly polygenic trait, given the many genes responsible for DNA repair, DNA replication, damage avoidance (e.g. anti-oxidant defenses), and mutagen detoxification, which in aggregate can determine MR. For example, many traits studied in humans have been shown to be highly polygenic [35]. The evolution of body size, somatic maintenance and germline mutation rate was then tracked under various regimens of selection (see also Methods: Model algorithm).
The model should approximate a sexually reproducing population. The model operates with single-parent reproduction model so that each individual descends from one parent. In this regard, technically it is tempting to view it as a model of an asexual population. However, the fundamental differences that lie between sexual and asexual populations (aside from the issue of purging deleterious mutations) are the amount of variation produced per the same sized population per generation and the presence of genetic recombination. Variance of inheritance in our model is too high to be assumed as having been generated by mutations accumulating along a clonal lineage and equals 10% of a trait's value per generation within 95 percentile. Genetic recombination associated with sexual reproduction segregates different genes and thus obstructs genetic linkage characteristic of inheritance during clonal reproduction. However, as Fig. 2 demonstrates, the role of genetic recombination in the evolution of multigenic traits is inversely proportional to the number of genes and alleles encoding a trait. The gene pool shown in Fig. 2a encodes a multigenic trait, whose distribution in a population is shown in Fig. 2c. As can be seen in Fig. 2c, an increasing number of genes encoding a trait decreases the variance of the encoded trait when randomly sampled from a population. This means that randomizing factors, such as genetic recombination, that randomly shuffle alleles for each gene in a population will have a smaller effect on the phenotypic expression of a trait encoded by many genes compared to oligo-and monogenic traits. Multigenic traits also provide multiple combinations of alleles for each of the encoding genes that can result in the same net expression of the phenotype. Based on this logic and simulation, mutator phenotypes, given the highly polygenic nature of mutation rate, will not be easily segregated from adaptive traits by genetic recombination as occurs with monogenic traits. We therefore did not simulate the processes of allelic segregation by recombination in order to reconstruct a sexual population (which would require many unwarranted assumptions), but modeled the distribution of mutation rate that was independently inherited from other traits and was not directly subjected to selection in the model (Fig. 2d). As such, the model only operates with the net ultimate change of a trait over generations. The assumed multigenic nature of the simulated traits also means that both segregation of alleles by recombination and aggregation of alleles by co-selection are impeded. We therefore assume that the net coevolution of a pair of multigenic traits will ultimately depend on selection acting on these traits that can overcome the effects of genetic recombination. Of course, in the real world, an individual could acquire a gene allele that substantially increases gMR, and we would expect that purifying selection would limit the ability of this allele to persist in the population (given its fitness cost); however, we are modeling under the assumption that the effect of most mutations will have an incremental effect on MR. We acknowledge that our modeling will be limited by our assumptions for how MR is encoded. While we explicitly model the cost of sMR realized in reduced body fitness, our model does not incorporate lethal germline mutations, as they affect the population randomly (all zygotes have roughly the same chance of acquiring a lethal mutation) and thus are irrelevant to selection. Our model operates with distributions of phenotypic effects, DPhE, of germline mutations, which in general can be safely modeled as normally distributed in accordance with the Shelford law of tolerance [36] which defines a typically normally distributed fitness of a phenotype within the tolerable range of environmental factors. The cost of germline mutations (which relates to At each timepoint during the simulation the modeled population undergoes 5 main updates: 1. Individuals that have not reached maturity increase their body mass following their growth curve, starting from the initial birth mass and up until they reach their inherited body mass (parent body mass with variation proportional to parent mutation rate); mature individuals remain at the same body mass. 2. Each individual past maturation reproduces with a certain inherited frequency of reproduction, producing on average an inherited number of progeny per litter, each progeny's birth body mass is inherited from its parent with variation proportional to the parent's mutation rate. Each individual is tried in a binomial trial with a small probability (at each timepoint) of dying of three main causes: 3. Death following limitations imposed by ecosystem carrying capacity which allows for a certain maximum population size and promotes intra-specific and inter-specific competition for resources if population numbers exceed this capacity. 4. Death caused by predation is modelled based on the Lotka-Volterra model of predator-prey interaction [34]. 5. Death caused by physiological aging, such as due to cancer, frailty or other age-related causes; the probability is negligible early in life but increases exponentially with age; the speed of increase of the probability of death caused by aging depends on an individual's aging profile which is determined by the aging curve as explained in Fig. 3a, "Theoretical introduction to the modeling" subsection of Results and "The somatic maintenance program paradigm" subsection of Methods the distribution of fitness effects, DFE) depends on and is imposed by the mode of selection applied to a trait as shown in Fig. 2e. For instance, under stabilizing selection (population mean trait is favored), a high germline mutation rate which produces more widely distributed trait values (wider DPhE) will impose a fitness cost on mutations with large phenotypic effects (deviating from the mean). Under positive selection, respectively, the positively selected tail of the DPhE (deviant trait values) will be favored. However, under both regimens the vast majority of phenotypically expressed germline mutations will lower fitness and be in the negative tail of the DFE (Fig. 2e). Our model for the cost of germline mutations therefore avoids artificial assumptions and replicates the . c A Matlabgenerated (see code in Additional file 1 Section 6) distribution of phenotype expression in a population of 100,000 individuals with different number of genes encoding the trait and each represented by 8 alleles. The X-axis in all the three charts is the same from 1 to 8, the Y-axis is population frequency of the phenotype (no scale for simplicity). Arrows demonstrate the allelic composition of phenotype examples generated based on multigenic (panel a) and monogenic (panel b) trait inheritance. The probability that all alleles in a multigenic trait will change the trait in the same direction (for instance all A8 alleles in panel a) is inversely proportional to the number of genes and is the product of the probabilities of each particular allele occurring in the given phenotype. An increased number of genes affecting the trait, therefore, decreases the likelihood that a random sampling event, such as genetic recombination, will result in a significant deviation of the resulting trait expression. d Assuming the multigenic nature of mutation rates and body size (primary modeled traits), we modeled both as normally distributed (see panel c) and independently inherited. Selection was applied on body size (or other modeled trait) but not on mutation rate. The evolution of mutation rate in the model was completely independent of other traits. e The left and right charts demonstrate a typical distribution of phenotypic trait expression (DPhE) in a population. Directional selection (left chart, green arrows) will favor extreme phenotypes from the target phenotype tail (left chart, green box). Stabilizing selection (right chart, green arrows) will favor the population mean which is the target phenotype (right chart, green box). Under both scenarios, most germline mutations will be in the negative (disadvantageous) tail of the distribution of fitness effects, DFE (middle chart), imposing a fitness cost for germline mutations. The exact ratio of the tails of the resulting DFEs (and thus the cost/benefit ratio of germline mutations) are determined by the mode and strength of selection acting on a phenotype natural mechanism for how selection determines the fitness cost of germline mutations. Noteworthy, this cost will differ for different traits depending on the current DPhE of the trait in the population at any given moment of evolution, as well as the current mode and strength of selection acting on the trait.
Neither equivalence of nor tight correlation between sMR and gMR are assumed in the model. We assume that sMR and gMR are linked by the underlying shared cellular machinery. However, this linkage is not expected to be rigid in nature. The model's "mutation rate" parameter does not produce equal sMR and gMR. The effective sMR in the model is generated by assigning the "mutation rate" parameter a fitness cost, which impacts the amount of actual damage to the soma resulting from mutations and thus acts as a proxy for the natural sMR, and the gMR is generated by assigning the "mutation rate" parameter a value of its effect on variance in the progeny. Therefore, the "mutation rate" parameter in the model is an overarching interface (the linkage factor) that is then independently realized as sMR and gMR in their respective tuning parameters (somatic fitness cost and progeny variation, respectively). We would add that mutations refer not just to single base changes (as often quantified when mutation rate is described), but all heritable changes in the DNA code, including deletions, insertions, and chromosomal rearrangements. The rates of these larger changes across individuals and species are poorly understood, but our sMR and gMR should incorporate all such heritable changes.
The model incorporates three major factors of mortality, including aging. Human life tables indicate that aging proceeds exponentially, whereby mortality and diseases accelerate at advanced ages (e.g. https://www.ssa.gov, https://seer.cancer.gov). The combined action of SMP mechanisms provides for an extended early period of high body fitness with minimal decline. We generalized this complex program in a curve that describes modeled animal mortality of physiological causes schematically shown in Fig. 3a and based on the following equation: where D A is the probability of dying of physiological causes at age A, M is mutation rate, and Som is a composite parameter that determines SMP efficiency. The cumulative distribution function of D A , or the probability of dying of physiological causes by age A, resembles human mortality (Fig. 3b). Note that physiological causes of death could include multiple aging-associated disorders, including cancer. The equation should thus provide a robust model for aging-related mortality, reflecting the extended period of high fitness and the late-life accelerating mortality. Figure 3a also demonstrates the relative effects of MR, which is a linear contributor, and the Som parameter, which stands for the total damage buffering capacity of the SMP (for details and theory see Methods: The somatic maintenance program paradigm). It is important to keep in mind that the M parameter (mutation rate) in Eq. 1 is responsible for the somatic costs of MR (higher MR in Fig. 3a accelerates aging-related mortality). Improved SMP, just as body size or other traits, may come at a cost on a short evolutionary time scale, which is later diminished by further adaptation. We did not include this cost in the modeling, since if a trait responds to directional selection this means that the benefit outweighs the cost. Since the amount of change of a trait as a result of positive selection in our model is arbitrary, we can conclude that this amount of change could already incorporate the net benefit minus cost. In other words, if the benefit of an evolutionary change exceeds its cost, then modeling benefit and cost on an arbitrary scale is mathematically equivalent to modeling only benefit. A summary of life history parameters can be found in Table 1. Thus, individuals in the simulated population inherit a number of life history traits, such as body size, longevity, body mass at birth, litter size etc. Each individual also inherits the mutation rate. Each of the inherited traits is inherited with variation proportional to mutation rate. For example, the inherited parental value for adult body mass is X. Let V be standard deviation of inheritance that is proportional to mutation rate. Then the progeny will grow to the size (X + V) drawn as a random number from a normal distribution with the mean X and standard deviation V. In other words, the higher the mutation rate, the more variant progeny's traits become relative to those of their parents. The mutation rate trait is also itself inherited with variance that is proportional to the current mutation rate (to itself; variance proportional to the mean). Therefore, the higher the inherited parental mutation rate, the more variance demonstrates the inheritance of mutation rate by progeny.

The evolution of SMP promotes selection for increased body size through better tolerance of MR
In our simulations, positive selection for increased body size (Fig. 3c, green) led to a concurrent selection for elevated gMR (Fig. 3d, green) and improved SMP (Fig. 3e, green). Artificially blocking SMP evolution by fixing SMP at the initial value (Fig. 3e, blue) significantly slowed the evolution of body size (Fig. 3c, blue; p < < 0.001) and triggered selection for lower gMR (Fig. 3d, blue). We implemented the ecosystem carrying capacity by setting a maximum biomass for the population; therefore, increasing body size led to a corresponding decline in population numbers, amplifying the power of drift ( Fig. 3f, g). When SMP was allowed to evolve, however,  the population entered a "drift zone" when its size decreased to~4000 individuals, which shortly thereafter was overcome by selection for even larger body size, visible also by a continuing decline in population numbers (Fig. 3f). When we artificially blocked the evolution of SMP, however, the drift zone was more profound. It occurred earlier at the population size of~6000-7000 individuals, and the population was not able to escape from it (for~1000 generations) and restore its initial rates of evolution (Fig. 3g), indicating an important role of SMP evolution in maintaining evolvability. We further generated a population with two simulated genotypes -Genotype A that could evolve SMP (10% of the population) and Genotype B with SMP fixed at the initial value (90% of the population). We set a maximum population size and removed the maximum biomass limit to rule out body mass effects on population size and selection, and tracked Genotype A and Genotype B frequencies under positive selection for increased body size (for code see Additional file 1: Section 1b). Despite the initial abundance, Genotype B (with fixed SMP) lost the competition in less than 200 generations, reflecting a direct competitive advantage of the capacity to evolve enhanced SMP (Fig. 3h). Hereafter, we will call the setting with positive selection for increased body size and freely evolving SMP and gMR the standard condition (usually shown in green, unless otherwise indicated) used in comparisons with other selection regimens.
Abrogating selection for increased body size reduces selection for gMR and SMP In the absence of positive selection for increased body mass (Fig. 4a, blue), both gMR (Fig. 4b, blue) and SMP (Fig. 4c, blue) demonstrate early positive selection, which appeared to have been caused by rapid evolution of reproductive parameters (see Additional file 1: Section 2). Overall, gMR demonstrates a significant general decrease [non-overlapping confidence intervals (CIs) at the beginning relative to the end of the simulation], and SMP undergoes a significantly smaller improvement compared to the standard condition (green; p < < 0.001). Blocking the evolution of body mass (Fig. 4d, blue) and SMP (Fig. 4f, blue) expectedly led to the evolution of lower gMR (Fig. 4e, blue) compared to the standard condition (p < < 0.001), which we interpret as being driven by the sMR costs in the absence of benefits of high gMR. In other words, mutation rate is selected against because of its somatic costs and the absence of benefits of higher gMR in static conditions. In natural populations that are under stabilizing selection, gMR will have costs due to greater phenotypic variance from a well-adapted state that are independent of sMR.

Decoupling sMR cost from gMR enhances the evolution of larger bodies
To investigate the role of the putative gMR benefit versus sMR cost balance in evolution, we further decoupled gMR and sMR by allowing gMR to evolve but making sMR cost fixed and independent of gMR (see Methods: Model variations). Decoupling sMR cost from gMR significantly accelerated the evolution of body size (Fig. 4g, blue) relative to the standard condition (green; p = 0.0052), revealing that sMR costs can limit the evolution of larger body size. During the early fast evolution of body mass, gMR (Fig. 4h, blue) and SMP (Fig. 4i, blue) demonstrate a corresponding positive response. Later, further body mass evolution becomes impeded (likely because of the severe depletion in population numbers), coinciding with the evolution of lower gMR. SMP plateaus during this second phase at a significantly lower level compared to the standard condition (p < < 0.001), indicating that the somatic costs of mutation rate stimulate the evolution of more robust SMP.
Selection acting on a somatic cost-unrelated trait still promotes the evolution of increased gMR and SMP As we have seen under blocked selection for increased body size (Fig. 4b, c, blue), SMP demonstrates an early phase of positive selection (Fig. 4c, blue) that is apparently reflected in a corresponding positive selection for higher gMR (Fig. 4b, blue). This observation suggests that both SMP and gMR may also respond to selection acting on some other traits, e.g. reproductive parameters (Additional file 1: Section 2). This raises the question whether SMP and gMR evolution would be sensitive to strong selection for a trait that does not affect somatic risks (greater body size increases the target size for somatic mutations). We simulated a condition that was similar to the standard condition, except positive selection was applied to a trait that did not affect sMR related somatic costs (see Methods: Model variations); e.g. if SMP improvement is solely a response to the increased sMR cost imposed by larger body, selection acting on an sMR cost-unrelated trait should not drive improvements in SMP. As shown in Fig. 4j (blue), unimpeded by increased sMR costs and declining population size, the evolution of an sMR cost unrelated trait is significantly faster compared to the evolution of increased body size (p < < 0.001). Interestingly, gMR (Fig. 4k, blue) also demonstrated an early phase of positive selection during early rapid evolution of the selected trait and remains above the initial gMR throughout the entire simulation. As anticipated, SMP is positively selected; however, in the absence of an increasing sMR cost associated with larger bodies, SMP's improvement is significantly smaller (Fig. 4l, blue, p < < 0.001). Notably, even with much less enhanced SMP, gMR is still under positive selection in response to positive selection for the sMR cost unrelated trait (Fig. 4k, blue), consistent with the sMR/gMR cost/ benefit ratio being an important factor regulating selection acting on gMR. Regardless, the results demonstrate that both gMR and SMP are responsive to selection for somatic risk unrelated traits, which indicates that high mutation rate is beneficial in positively selective conditions.

SMP enables maintenance of gMR when directional selection is weak
As we have seen in Fig. 4d-f, in the absence of strong positive selection for increased body size and SMP efficiency, selection acts to lower gMR. Figure 5 shows, however, that this selection is significantly modified by the efficiency of SMP. Stronger SMPs (lower Som value) relax selection for lower gMR when directional selection is weak (non-overlapping CIs between the standard (red) and either of the improved SMPs). If SMP is strong enough, early evolution for higher mutation rate is observed (Fig. 5, green and blue curves); as we explained this occurs because other traits such as reproductive parameters are evolving during the simulation for some period of time until they reach stasis (as previously shown in Additional file 1: Figure S1). After the stasis was reached, lower gMR is selected for independent of the strength of SMP, as it no longer provides any selective benefit. As will be explained further below, this observation may have significant implication on longterm species survival in relatively static environments where stabilizing selection dominates.

Modeling competition between wildtype and mutator phenotypes
Under strong positive selection, whether for increased body mass (Fig. 3a-c, blue) or a sMR cost unrelated trait (Fig. 4h, i, blue, and Fig. 4k, l, blue), we observed consistent signs of positive selection for higher gMR. However, because gMR and sMR are linked, and gMR also confers a cost at the organismal level by increasing phenotypic variance, higher gMR is a trait that should negatively impact individual fitness and therefore be under negative selection. To investigate this question, we mixed two simulated genotypes, one "wild-type" (50%) and one "mutator" (50%) in a population of stable size and under positive selection for a sMR cost unrelated trait. We then observed the genotypes' frequencies in the population using varying strength of mutators. Figure 6a demonstrates that while the mutator's fitness initially is lower compared to wild-type, eventually the mutator outcompetes its wild-type counterpart. Interestingly, with increased mutation rate, the magnitude of the mutator's initial decline increases, but so does the speed at which it subsequently overtakes the population. This result provides a clue for how higher mutation rate, being a trait with negative impact on fitness, can be selected for. Because net organismal fitness is a composite trait impacted by the fitness value of many individual traits, the initial fitness of the "mutator" is lower because, all other traits equal, higher MR incurs increased sMR cost. However, in response to selection, a mutator is capable of more rapidly developing other (adaptive) traits (Fig. 6b) and thus its overall fitness soon becomes higher compared to wild-type. Noteworthy, in Fig. 6 we modeled just two genotypes (wildtype and mutator) in competition starting from arbitrary initial phenotype frequencies so that the initially disadvantageous mutator phenotype has a high enough frequency not to be lost stochastically, and we do not model low-frequency phenotypes. In natural populations, mutation rate will represent a distributed trait encoded by many genes and alleles so that selection will act on multiple phenotypes with various degrees of expression of higher mutation rate, and with a corresponding range of frequencies depending on where a particular phenotype is in the distribution. In other words, there will always be phenotypes that have a sufficiently high frequency to prevent stochastic loss. In a monogenic trait setting, however, with two distinct phenotypes, the initial low frequency of a mutator will pose a risk of stochastic loss, which is beyond the subject of our study.

Discussion
Our study demonstrates that positive selection for increased body size triggers a concurrent selection for improved somatic maintenance to mitigate the increased somatic risks of larger bodies. Improved somatic maintenance, in turn, promotes the evolution of higher germline mutation rates by reducing the cost of somatic mutations and thus altering the sMR/gMR cost/benefit ratio. Conditions of strong positive selection for somatic cost independent traits, as our model shows, can also ; a linear decrease in the Som value results in a substantially improved SMP, so that the green SMP is~10X more efficient compared to red, and the blue is a~4X more efficient SMP than the green. The standard (red) SMP leads to a significantly stronger selection for lower gMR (non-overlapping 95% CIs); however, the absence of difference between the 10X (green) and 40X (blue) improved SMPs indicates that overly improved SMPs might not provide any further difference for how selection acts on gMR. Note also that even with enhanced SMP (green and blue) there is significant decrease (non-overlapping 95% CIs) in gMR relative to either the initial gMR or to the gMR at the peak (resulting from selection to optimize other traits) alter this balance by elevating the benefits of higher gMR. Alternatively, under stable conditions the sMR/ gMR cost/benefit balance is altered by the existing cost of somatic mutations and by the increased cost and absent/reduced benefits of higher gMR itself, which ultimately favors lower mutations rates. Under stasis (stabilizing selection dominates), gMR exerts a cost independent of somatic risks by increasing deviation of progeny phenotypes from population mean/median and thus reducing their fitness. Our study thus demonstrates that the evolution of mutation rate is not under a universal population size-dependent selection acting to lower it, but is highly tunable and governed by selection acting on other traits. Importantly, our modeling indicates that under certain conditions elevated mutation rate, unlike perhaps any other trait, can be positively selected despite its negative effects on individual fitness, as a trait (like body size) directly under positive selection is more likely to arise in an individual with higher MR. Mutation rate in eukaryotes is a highly polygenic trait encoded by multiple genes involved in DNA replication, repair, damage avoidance, and cell division machineries (see Introduction). Animals mostly reproduce sexually, which should generate an extensive population allelic diversity for these genes. This diversity should provide for a relatively continuous distribution of mutation rate in populations, rather than being a uniform trait marked with sporadic monogenic mutants, as may occur in asexual populations [37][38][39]. Such intra-population variation [40,41], as well as the ability of mutation rate to rapidly evolve [42] (at least for particular trinucleotide contexts) has been shown for humans. What drives or maintains this variability (if "driven" at all), whether in humans or other animals, is currently not established. Our modeling would suggest one possible contributing factorthat different populations have experienced differences in directional selection in the relatively recent past. Of course, there are other potential explanations, including bottlenecks and geographic isolation.
Sexual reproduction would be supposed to effectively segregate alleles contributing to mutation rate from alleles for other (e.g. adaptive) traits. It has been argued based on other evidence that the efficiency of such segregation in sexual populations is limited [43]. In particular, as we have argued in Theoretical introduction to the modeling and shown in Fig. 2a-c, that the multigenic nature of the gMR trait should substantially reduce the effects of genetic recombination on the evolution of the net gMR phenotype and its segregation from other traits. For any given allele that has a minimal impact increasing gMR, selection will have a reduced ability to eliminate this allele from the population. But selection could still act, albeit slowly, on the composite MR phenotype encoded by a large number of genetic variants.
It also appears from our results that animal evolution, with the macroscopic trend toward larger bodies, should have driven a concurrent evolution of extended longevity, the latter being determined by the efficiency of species-specific somatic maintenance programs. Even though extended longevity tentatively appears to be a benefit on its own, e.g. due to extended reproduction period, our model demonstrates that somatic maintenance (and thus longevity) is under a much weaker positive selection in the absence of other positively selected traits. This observation can explain why extended longevity demonstrates significant deviations across animal taxa from the general rule larger body → longer lifespan. Our results indicate that the evolution of longevity (as a function of somatic maintenance efficiency) should be greatly impacted by the rate of evolution of other traits, and not necessarily body size. Of course, longevity will also be highly regulated by external hazards and other extrinsic factors influencing survival probabilities, which has substantial experimental and modeling support [44], and which were not explored in our study.
Interestingly, our study predicts an important evolutionary role for the mechanisms of somatic maintenance in addition to their evolution as a means of improving individual survival of large animals [25,32]. Our results demonstrate that selection for enhanced somatic maintenance goes well beyond the evolution of body size and is promoted by strong directional selection acting on any trait. This result indicates that SMPs may have had an important role in the evolution of large animals. Selection for higher gMR following improved SMP may be an important mechanism "rescuing" the reduced responsiveness to selection imposed by reduced population size, extended generation times and lower reproduction rates in large animals. Therefore, SMPs and longevity may have an important contribution to species' long-term survival. For example, prolonged periods of relative evolutionary stasis [45][46][47][48] should trigger selection for lower mutation rates. By relaxing selection for lower mutation rate and thus maintaining evolvability (as shown in Fig. 5), enhanced SMPs can ensure better survival of animal groups facing rapid evolutionary transitions or drastically changed environments after such relatively static periods. All other traits equal, species with extended longevity may survive such transitions with higher probabilities.
Lynch and colleagues have provided extensive arguments supporting the idea that the higher MRs in animals compared to unicellular organisms are likely to be caused by reduced population sizes that limit the ability of selection to lower mutation rate [10][11][12]. In conjunction with population size, in large animals the strength of selection will be further attenuated by lower reproduction rates and extended generation times. Based on our results, Lynch's theory can be extended by recognizing that somatic maintenance programs (and longevity) should have substantial influence on the general relationship between population size and mutation rates, and on the strength and directionality of selection acting on mutation rates. For example, in our simulation, populations of the same initial size but with different SMP efficiencies demonstrate profound differences in the effects of population size-driven weakening of selection (Fig. 3f, g), as well as discrepant selection for mutation rates (Fig. 3d).
It is unlikely that DNA maintenance is "as good as it gets" [49]. Mutation rates vary widely across organisms [10], suggesting that mutation rate is not as low as theoretically possible and may be subject to selection and/or drift. Mutation rate could theoretically approach some lower limit under conditions where stabilizing selection dominates; but this situation could change during periods of directional selection (such as following environmental change), and our modeling indicates how the evolution of SMP could influence the evolution of the selected traits and mutation rate. Regardless, the reality is that we simply do not fully understand the limits of mutation rate. Even if standing variation may be sufficient under many circumstances to enable directional selection [50], our modeling suggests that baseline mutation rates could be higher than they otherwise would be given the evolution of SMP to mitigate their impact on the soma. Moreover, as shown in Fig. 7, selection for extreme trait values even from standing variation will likely co-select for individuals with higher mutation rates. Selection for higher mutation rates has been shown experimentally in bacteria [37][38][39]51], whereby engineered or spontaneous mutants with higher mutation rates have been shown to have advantages over wild-type under positively selective conditions. Mathematical modeling has been used to argue that stress responseregulated mutagenesis can facilitate adaptation to new environments [52][53][54]. The "mutator hitchhiker hypothesis" explains such selection by the higher probability that adaptive mutations will appear in a mutator cell [39]. Once such a mutation occurs, the mutator genotype spreads to fixation by being genetically linked to the adaptive phenotype.
Modeling studies demonstrate that evolution of evolvability, including varying selection on mutation rates, should be possible in sexually reproducing organisms [43,55,56]. Yet robust experimental corroboration of such a possibility appears to be lacking. Higher mutation rate exerts negative fitness effects on the current generation in any species by imposing somatic risks and contributing to aging. However, in changing conditions whereby directional selection prevails, higher germline mutation rates should provide for better survival of progeny, since directional selection favors phenotypic extremes (tails) in the distribution of phenotypes in a population (Fig. 7). Such extreme phenotypes are likely to descend from parents with higher mutation rates (relative to population mean) and thus inherit enriched pools of mutator alleles for multiple MR-related genes. While genetic recombination should be effective in separating any given allele from the adaptive alleles (those that are actually selected for), it should be less effective in segregating the entire multigenic mutation rate trait from the adaptive phenotype (see Fig. 2a-c). Under stabilizing selection whereby population mean (non-deviant phenotypes) is favored, a similar mechanism should trigger co-selection for lower mutation rate by repeatedly sampling the least deviating pool of phenotypes in the population. For conditions favoring either directional or A B Fig. 7 A theoretical model of the effect of selection regimens on the evolution of mutation rate. a Under conditions of stabilizing selection, phenotypes less deviant from the population mean in their adaptive traits (left Y-axis for frequency distribution) are favored; such phenotypes should inherit and pass on over generations lower mutation rates (sampled from the mutation rate distribution of the population) providing for smaller variation among progeny and thus better survival. Stabilizing selection leads to co-selection for phenotypes closer to the population mean and lower mutation rates. Improved somatic maintenance, by lowering the costs of MR, permits greater persistence of phenotypic variability during stasis. b In the conditions of major evolutionary transitions whereby positive selection dominates, the favored phenotypes for adaptive traits are in the selected tail of the population's phenotypic distribution; such phenotypes should inherit and pass on over generations higher mutation rates from the mutation rate distribution of the population. Positive selection leads to co-selection of phenotypes that deviate from the population mean and higher mutation rates. Improved somatic maintenance, by lowering the cost of MR, will enhance selection for altered trait values that will have increased prevalence in individuals inheriting higher MR (farther right on the schemata) stabilizing selection, investment in somatic maintenance will reduce the somatic costs of MR and thus increase selection for or tolerance of (respectively) MR.

Conclusion
Our results raise the question of whether the evolution of large body size in animals would be possible without such a complex pattern of selection acting on mutation rate, and whether such a complex relationship is necessary to explain the evolution of large animals. The evolution of large bodies has entailed the cost of losing the ability to evolve via all major parameters that define this ability, such as population size, reproduction rate and generation time, except mutation rate (which increased). Therefore, one scenario could have been that this cost has been so prohibitive for many species that positive selection for mutation rate was necessary to allow evolution of large animals. Alternatively, mutation rate could have been high enough to maintain evolvability at the selection/drift barrier point where selection was no longer able to reduce it further [10]. Understanding which of these scenarios prevails in the evolution of large animals requires more research.

Software
The model was created and all simulations were run in the Matlab environment (MathWorks Inc., MA) version R2014a.

Model algorithm
The model is a stochastic Monte Carlo type model (the exact algorithm can be found in Additional file 1: Section 1a) that runs a total of 1,005,000 updates ("time" in arbitrary units, AU) unless otherwise stated, which represents~1000 generations of the simulated animal population (see Fig. 1 for the flow chart). The simulation starts with building an initial population of 10,000 individuals. Each individual has a number of simulated traits: 1) ID, which is 1 (monogenotypic population) or 1 and 2 (in experiments with competition between two genotypes in a mixed population to indicate genotypes); 2) current age, which increments by 1 at each simulation update; 3) inherited body mass, which is inherited with variation by an individual and will be reached by adulthood (at age~1000) and equals 5000 AU in the initial population; 4) current body mass, which changes during individual growth, following a growth curve, and plateaus at the inherited body mass in adults; 5) inherited birth mass, which in individuals of the initial population is 300 AU; 6) inherited mutation rate of 10 − 9 AU (explained below); 7) inherited reproduction rate, which is the period with variation between successive reproductions in adult individuals and equals~600 in the initial population; 8) inherited litter size (initially 1), which is the number of progeny produced per individual per reproduction; 9) inherited parameter of somatic maintenance, which determines the strength of the somatic maintenance program as further explained below; 10) age of first reproduction, which dictates that an individual begins reproducing when its current body mass reaches 0.9693 of its inherited adult body mass (the number is derived so that in the initial population maturity is reached at age~1000 based on the growth curve).
Each inherited trait varies in progeny relative to parental. This variation was produced by multiplying the inherited mutation rate by the parameter of inherited variance (inhvar = 250,000,000) and the product was used as the standard deviation (STD) of the normally distributed variation in inheritance. This transformation was not necessary, as the inhvar parameter is constant throughout simulation and it simply determines the magnitude of the mutation rate's effects in the germline, which is imaginary and in the initial population simply produces 0.000000001 × 25,000,000 = 0.025 that serves as the STD parameter for the normal distribution from which inheritance variation is drawn. However, we kept this two-parametric model for inheritance because mutation rate is also separately used in the equation of the somatic maintenance program (as will be explained later).
Each newborn individual grows, reaches maturity, then reproduces over the rest of its lifetime and eventually dies. The model is asynchronous, so that at every time-point of the simulation the population contains individuals of various ages whose lifecycles develop independently.
And finally, three factors of mortality were modelled in the simulations. First, at every timepoint of the simulation, an individual could die of somatic causes with a certain probability. This probability is small at the beginning of life (but still can be caused by some imaginary inherited genetic defects) and increases exponentially with age based on the paradigm of the aging curve, which is primarily determined by an individual's inherited somatic maintenance program (SMP). In humans, the aging curve also depends on lifestyle, however we assume in this model that in a wild animal population lifestyle distribution is sufficiently uniform to be neglected. More detailed description of the somatic maintenance paradigm that we applied will be explained further below. Secondly, the simulated animals had a chance of dying of external hazards, such as predators. We applied the Lotka-Volterra model of predator-prey interactions [34,57] to implement the dynamics of predator pressure (effectively the chance of dying of an external hazard cause per time unit). Here we should mention that smaller individuals and juveniles had higher chances of dying of external hazards, which effectively created positive selection for increased body size and also reflected the typical high mortality rates among juveniles observed in natural populations. And lastly, individuals could die of intra-specific competition. We implemented such competition by setting the upper limit of population's total biomass, which in nature is imposed by the ecosystem's carrying capacity. Therefore, in the simulated population biomass produced over the biomass limit caused additional mortality, so that stochastically, population total biomass never exceeded the limit. Larger individuals also had lower probability of dying of intra-specific competition, based on the assumption that competition for resources and mates (the failure to reproduce is effectively an evolutionary death) will typically favor larger individuals and this should have been one of the forces that has been driving the macroscopic animal evolutionary trend towards increasing body size. The advantage of size in this mortality model also created additional positive selective pressure for body size. The total age-dependent mortality of all causes in our model did approximate a typical wild animal mortality curve (Additional file 1: Section 3).

The somatic maintenance program paradigm
In order to replicate natural mortality caused by physiological aging, such as cancer, decreased immune defense and lower ability to avoid predators or to succeed in intraspecific competition, we made use of the aging curve, or somatic maintenance, concept. Modern human (in developed nations) and captive animal mortality curves (Fig. 3b for human) differ from wild animal mortality curves in very high early life survival with most mortality significantly delayed into advanced ages [58,59]. This difference is caused by many reasons, such as much lower mortality caused by external hazards and better nutrition and general healthcare. It therefore can be assumed that the human and captive animal mortality curves are close representations of the physiological aging curve. As longevity depends on multiple mechanisms of maintaining the soma, we can also call this curve the somatic maintenance curve. In order to reconstruct this curve, we assumed that somatic maintenance depends on the interaction of two opposing forces: 1) the accumulation of genetic and structural damage in the soma that promotes aging and 2) the somatic maintenance program consisting of a number of mechanisms that prevent or buffer the effects of genetic and structural damage. The exact mathematical relationship between these two forces and age is not known, however an example of cancer development can be used as a proxy to explain the equation we derived for it.
Oncogenic mutations (including oncogenic epigenetic changes) are the ultimate necessary condition for cancer to develop. The frequency of oncogenic mutations linearly depends on mutation rate on a per cell division basis. Therefore, we assume that linear changes in mutation rate will have linear effects on the odds of the occurrence of oncogenic mutations. An oncogenic mutation provides the initiated cells with a linear change in their fitness relative to normal cells. However, over time an advantageous clone with a constant linear fitness advantage will proliferate exponentially. Therefore, we can already assume that mutation rate should have a linear effect on the cancer curve, while time/age adds an exponential component revealed in an exponential growth of a tumor. We can reasonably assume further that a strong SMP will efficiently suppress such a clone, slowing or even preventing its growth [60]. A weaker SMP will allow the clone to proliferate faster. Therefore, SMP strength can modulate the effects of mutations and time on cancer risk. The exact relationship between SMP strength and physiological risk factors is not known. However, we know that their interaction leads to a net exponent in physiological decline and disease risk.
We therefore reconstructed the human aging curve by maintaining the general principal relationship between these factors as shown in Eq. 1. As seen from the equation, mutation rate is a linear contributor to aging. Age itself contributes exponentially, and the somatic maintenance composite parameter Som is, in turn, in power relationship with age. The cumulative distribution function of D A (Eq. 1) produces D(A)the probability of dying of somatic/physiological causes by age A and yields a shape close to the human mortality curve (Fig. 3a, b). We cannot claim that these three factors are in the exact relationship predicted by Eq. 1, as it is unknown. As seen in Fig. 3a, changes in the Som parameter have substantially greater effects on the resulting mortality curve than mutation rate, with mutation rate still having a sizeable effect as well. Yet claims are still made (e.g. [61]) that mutation rate is a larger factor in aging than we assume in this model. Validation of our assumption in general comes from the body of evidence that up to 50% of mutations in humans accumulate during body growth by the age 18-20 [62][63][64], although this early life pattern of mutation accumulation is not universally observed [65,66] (still, note that Osorio et al. quantify a~4-fold higher mutation rate in hematopoietic stem cells during fetal development relative to post-natal life). If mutation accumulation had a significant effect on aging on its own, we should age rapidly until age 18-20 (half-way) and then the rate of aging should decelerate. However, in reality the rate of aging (physiological decline) is highest after 50 years of age, indicating that the combined strength of the SMP has an overpowering effect in modulating the effects of genetic damage on aging. As a result, we reason that Eq. 1 might reasonably approximate the natural relationships among these three factors. Therefore, based on an individual's aging curve we calculated the D A parameter at each simulation time-point (using the individual's mutation rate, age and Som parameter) and applied it in a binomial trial as the probability of that individual's dying of somatic/physiological causes in an age-dependent manner. As further explained in Additional file 1: Section 4, the exact relationship between the Som parameters and each of the other two (mutation rate and age) has no effect on the model, as the model represents SMP and its variation by using area under the mortality curve. Therefore, the sole purpose of Eq. 1 in the model is to generate an age-dependent curve of physiological mortality whose cumulative function (probability of dying by a certain age) resembles in shape the human mortality/ aging curve (see Additional file 1: Section 4 for detailed explanation and illustration).

Model variations
A number of model variations used in simulation experiments are employed. Fixed trait values involved simply fixing the initial trait value without inherited variation throughout the entire simulation. Dislinking of somatic and germline mutation rate was done by making the value M in Eq. 1 independent of an individual's mutation rate, which resulted in somatic costs becoming independent of transgenerational variation of mutation rate (effectively from germline mutation rate). Selection for a trait that did not affect somatic risks was achieved by transforming the "body mass" trait's effects by removing the trait from calculations of the risk of death by somatic causes (unlike body size, it did not influence the risk), then removing the population biomass limit and setting maximum population size (unlike body mass, other traits do not directly affect population numbers) and fixing the growth rate curve so that it reached the initial body mass of 5000 AU (the current body mass parameter in the model; the inherited body mass variation did not exist and the inherited body mass parameter was replaced with the somatic risk unrelated trait). These manipulations made the selected trait a proxy for a trait unrelated to somatic risks (e.g. hair color). Competitive assays included individuals with different ID parameters, such as 1 and 2 to indicate different "genotypes"; traits of the "genotypes" then were tracked and stored separately.

Data processing
Processing of primary data included removal of outliers (see Additional file 1: Section 5). Occasionally the simulations generated "NaN" (not a number) values in individual parameters, which were rare but quickly propagated if left in the population. We immediately deleted individuals from the population if "NaN" values appeared in any of their parameters. Based on the rarity of such events, we can assume that they had the effect of rare early lethal mutations and affected the population at random. Thus, we assume these did not affect the principal results.

Statistics and data presentation
Most simulation experiments were made with 25 repeats. Due to heavy skews in sample distributions (inferred by D'Agostino-Pearson test for normality of a distribution), all figure panels represent medians (thick lines) and 95 percentiles on each tail (color-shaded areas). Statistical differences between experimental conditions were calculated as follows. We first calculated the sum of all values in each run throughout the entire evolution of a trait (typically 1,005,000 time points). In this way, given the small increment over a long time the sum essentially approximated the area under the curve of a trait's evolution. These sums (usually 25 repeats in one experiment/sample) were then compared by applying the Matlab implementation of the Wilcoxon rank sum test, which is considered equivalent to the Mann-Whitney U-test. P-values <= 0.05 were considered as indicating significant difference.