Horizontal gene transfer (HGT) is a major factor in the evolution of prokaryotes. An intriguing question is whether HGT is maintained during evolution of prokaryotes owing to its adaptive value or is a byproduct of selection driven by other factors such as consumption of extracellular DNA (eDNA) as a nutrient. One hypothesis posits that HGT can restore genes inactivated by mutations and thereby prevent stochastic, irreversible deterioration of genomes in finite populations known as Muller’s ratchet. To examine this hypothesis, we developed a population genetic model of prokaryotes undergoing HGT via homologous recombination. Analysis of this model indicates that HGT can prevent the operation of Muller’s ratchet even when the source of transferred genes is eDNA that comes from dead cells and on average carries more deleterious mutations than the DNA of recipient live cells. Moreover, if HGT is sufficiently frequent and eDNA diffusion sufficiently rapid, a subdivided population is shown to be more resistant to Muller’s ratchet than an undivided population of an equal overall size. Thus, to maintain genomic information in the face of Muller’s ratchet, it is more advantageous to partition individuals into multiple subpopulations and let them “cross-reference” each other’s genetic information through HGT than to collect all individuals in one population and thereby maximize the efficacy of natural selection. Taken together, the results suggest that HGT could be an important condition for the long-term maintenance of genomic information in prokaryotes through the prevention of Muller’s ratchet.
Evolutionary advantages of recombination have been extensively considered, primarily in relation to the evolution of sex in eukaryotes (Kondrashov 1993; Maynard Smith 1978, 1998). However, advances in genome sequence analyses have revealed that many prokaryotes, whose mode of reproduction is asexual (i.e., no meiosis or syngamy), frequently undergo genetic recombination, the process known as horizontal gene transfer (HGT; Feil et al. 2000, 2001; Ochman et al. 2000; Koonin et al. 2001; Gogarten et al. 2002; Vos and Didelot 2009; Yahara et al. 2012). Numerous studies indicate that HGT has been instrumental in producing functionally consequential genomic changes throughout the evolution of prokaryotes (Lerat et al. 2005; Dagan et al. 2008; Treangen and Rocha 2011). For example, HGT is implicated in the evolution of various functional systems and specific adaptations such as various metabolic pathways (Pál et al. 2005; Maslov et al. 2009), oxygenic photosynthesis (Mulkidjanian et al. 2006), thermal resistance (Aravind et al. 1998; Brochier-Armanet and Forterre 2007), antibiotic resistance (Blahna et al. 2006), pathogenicity (Kelly et al. 2009), and others. Moreover, it has been shown that HGT can substantially impact microbial community development (Burke et al. 2011) and ecological diversification (Shapiro et al. 2012). The prevalence and significance of recombination in prokaryotes as well as in eukaryotes point to the crucial importance of recombination for the evolution of life in general (Komatsu 1980).
HGT in prokaryotes occurs though three basic mechanisms: conjugation, transduction, and transformation (Levin 1988; Thomas and Nielsen 2005; Johnsborg et al. 2007). Conjugation and transduction involve transfer of prokaryotic DNA between cells mediated by infectious agents such as viruses and conjugative plasmids, respectively. Because genes involved in these mechanisms reside in the genomes of mobile genetic elements, transduction and conjugation can be considered side effects of the selfish propagation of these infectious agents (Redfield 2001).
Transformation, by contrast, is a bacterium-programmed mechanism of HGT that is mediated by natural competence systems encoded in the genomes of many bacteria (Lorenz and Wackernagel 1994; Claverys et al. 2009). Through transformation, a cell absorbs and integrates extracellular DNA (eDNA) into its chromosome. Although transformation might be viewed as a side effect of metabolic consumption of eDNA (Redfield 2001), there is increasing evidence, at least in several bacteria, that absorbed eDNA serves as the source of information rather than or at least in addition to being a source of nutrients (Johnsborg et al. 2007; Claverys et al. 2009; Zafra et al. 2012; Johnston et al. 2013). Such evidence includes, e.g., the existence of competence-induced proteins (e.g., DprA in Streptococcus pneumoniae) that protect absorbed eDNA from degradation by intracellular nucleases (for other lines of evidence, see the references cited above).
Besides transformation, there is another mechanism of HGT called gene transfer agents (GTAs), which is encoded in the genomes of many prokaryotes (Lang and Beatty 2007; McDaniel et al. 2010; Lang et al. 2012). The GTAs are defective virus particles that encapsulate apparently random fragments of the host DNA (rather than the genes encoding the GTA proteins) and mediate the transfer of these fragments into other prokaryotic cells. At face value, the GTAs appear to function as dedicated vehicles for HGT. The existence of these apparently prokaryote-programmed mechanisms of HGT suggests an intriguing possibility that these mechanisms have been maintained by selection in prokaryotes owing to evolutionary advantages offered by HGT—a situation akin to what is currently believed for sexual recombination in eukaryotes (Michod and Levin 1988; Kondrashov 1993; Keightley and Otto 2006; see Bernstein et al. 1988 and Redfield 2001 for counteracting views).
What selective advantages could HGT offer to prokaryotes? To address this question, we asked whether an evolutionary advantage of recombination that has been originally suggested for eukaryotes could be extended to the case of HGT in prokaryotes. In particular, we hypothesized that HGT could restore genes that are deleted or inactivated by mutations, thereby preventing the stochastic, irreversible deterioration of genomes in finite populations known as Muller’s ratchet (Muller 1964; Felsenstein 1974; Pamilo et al. 1987; Bell 1988; Charlesworth et al. 1993; for other types of advantages, see, e.g., Redfield et al. 1997; Cooper 2007; Levin and Cornejo 2009; Vos 2009; Wylie et al. 2010). This hypothesis is consistent with the observation of extreme genome reduction and accelerated protein evolution in intracellular bacterial parasites and symbionts (Moran 1996; Novichkov et al. 2009; Merhej and Raoult 2011; McCutcheon and Moran 2012). Because of continual bottlenecks and sequestration within the host organisms, these bacteria are likely to have small effective population sizes and be effectively isolated from HGT between populations in different hosts. Under these conditions, the impact of Muller’s ratchet increases and could lead to genome reduction, accelerated protein evolution, and possibly eventual extinction (Lynch et al. 1993; Merhej and Raoult 2011). Therefore, HGT might be an important requirement for the long-term maintenance of genomic information in prokaryotes (Koonin 2011).
However, HGT in prokaryotes differs from eukaryotic recombination both in quantity and in quality. Quantitatively, HGT is relatively less frequent and affects a far smaller portion of chromosomes per recombination event than eukaryotic recombination that is inextricably linked to sexual reproduction (Dubnau 1999; Feil et al. 2000). Qualitatively, whereas eukaryotic recombination is a direct exchange of DNA sequences between haploid genomes inherited from the two parents, HGT in prokaryotes, in particular, transformation is an indirect exchange of DNA sequences between cells via an eDNA pool (Redfield 1988). The sources of eDNA can be diverse, but the most common one is likely to be the debris of dead conspecific cells. If the death rate of a cell increases with the number of deleterious mutations in the genome, the genomes of dead cells on average carry more deleterious mutations than those of living cells (Redfield et al. 1997). Consequently, transformation on average is expected to introduce more deleterious mutations than it would remove. In view of these differences between eukaryotic recombination and HGT in prokaryotes, which potentially nullify the hypothetical benefit of HGT, we sought to assess the net effect of HGT on the operation of Muller’s ratchet by performing stochastic population genetic simulations.
We analyzed a simple evolutionary model of finite prokaryotic populations undergoing HGT. Our model is based on the infinite population model studied by Redfield (1988) with several extensions (see the section Materials and Methods). The model was analyzed under the assumption that is most unfavorable to the hypothesized advantage of HGT, namely, that variation in fitness is attributed solely to variation in death rates, and the fittest individuals do not contribute their genomes to the eDNA pool. The results show that, even under this assumption, HGT can prevent the operation of Muller’s ratchet by causing a continual generation of the least-loaded class (i.e., genomes with the least number of deleterious mutations) from more-loaded classes. Moreover, the prevention of Muller’s ratchet by HGT is substantially more efficient in a subdivided population than in a panmictic population of an equal overall size, provided that the diffusion of eDNA is sufficiently rapid.
Materials and Methods
Definitions for the symbols used in this article can be found in Table 1.
General description of the model
The model consists of two components: a population of prokaryotic cells and a pool of eDNA. The population undergoes a standard mutation-selection process (Wright-Fisher model) as well as HGT, through which cells acquire alleles randomly drawn from the eDNA pool. The eDNA pool gains input from the population due to natural cell death and also undergoes spontaneous decay. The details of each of these processes are described in the sections to follow.
The population dynamics was modeled as the discrete-generation Wright-Fisher process with fixed population size N. The genome was modeled as l loci (haploid), each representing a genomic segment that could be replaced by one recombination event (Wylie et al. 2010). The typical lengths of such segments are known to be on the order of 10 kb in several bacteria (Dubnau 1999; Feil et al. 2000), and the lengths of bacterial genomes are on the order of 1 Mb. Accordingly, the value of l was set to 100 unless otherwise stated. Each locus was assigned one of the infinite number of possible alleles, each of which was represented by a pair of integers: one integer indicating the locus and the other one indicating the number of deleterious mutations. All deleterious mutations were assumed to cause an equal fitness effect s without epistasis; thus, the relative fitness fi of a genome with i mutations is fi = (1 − s)i. This simplified model was adopted so that the effect of HGT on Muller’s ratchet could be assessed without confounding the additional effects caused by epistasis (Redfield et al. 1997). The speed of Muller’s ratchet is known to depend on the two combinations of parameters, U/s and sN (Neher and Shraiman 2012). Thus, in this study, the value of s was set to 0.01, and the values of U and N (more specifically, sN exp(−U/s) (denoted as ) and N) were varied to explore the behavior of the model.
In the model, each genome acquired the Poisson distributed number of new deleterious mutations with the mean U per generation. Beneficial mutations were ignored because we were concerned with the possible deterioration or maintenance of initially well-adapted genomes, for which the frequency of beneficial mutations is likely to be rare (Maynard Smith 1978; Charlesworth et al. 1993; cf. Goyal et al. 2012).
Horizontal gene transfer
HGT generally can be categorized into homologous recombination and illegitimate (i.e., nonhomologous) recombination. Because homologous recombination is much more frequent than nonhomologous recombination (Hulter and Wackernagel 2008; Brigulla and Wackernagel 2010), only the former was included in the model. In each genome, homologous recombination replaced the Poisson distributed number of alleles with mean r per generation per genome with alleles randomly drawn from the eDNA pool (after the replacement the original alleles were discarded).
Dynamics of the eDNA pool
The frequency of homologous recombination decreases exponentially with the divergence between DNA sequences (Watt et al. 1985; Majewski and Cohan 1999; Eppley et al. 2007). Thus, the most important source of eDNA for homologous recombination would be the genomes of conspecific cells. For the sake of simplicity, the model assumed that genomes of conspecific cells in the same population were the only source of eDNA.
Prokaryotic cells can release their genomic DNA into the extracellular environment either actively by fratricide, suicide, or DNA secretion, or passively by natural death. However, the model assumed that genomic DNA was released only by natural cell death in order to consider the simplest and most unfavorable condition for the putative advantage of HGT (see the Introduction section).
To incorporate the influx of eDNA as the result of natural cell death, the model assumed that the per-generation probability of an individual with i mutations releasing its genome into the eDNA pool (i.e., dying) was where nj is the number of cells with j mutations, and is the fitness of cells with j mutations normalized by the fitness of the least-loaded class, i.e., where mLLC is the number of mutations in the least-loaded class (mLLC is defined as the integer for which nm = 0 for m < mLLC and ).
The aforementioned assumption on εi was obtained from the following heuristic argument. Let us consider a continuous population dynamics described by dni/dt = bni − dini, where b and di are the birth and death rates, respectively. The absolute fitness per unit time of a cell with i mutations can be expressed as . The relative fitness is expressed as to be consistent with the aforementioned definition of fi (so that f0 = 1). The fraction (probability) of cells dying per unit time is , which can be expressed as 1 − (1 − D0)fi, where (D0 is the probability that an individual of the least-loaded class dies). Thus, the normalized contribution to the eDNA pool by a cell with i mutations is . Setting D0 to zero maximizes the number of mutations in the genomes of dead cells (see Supporting Information, File S1), leading to the expression . This expression is used in the seminal work of Redfield (Redfield 1988), which considers an infinite population model. In a finite population model, however, the number of mutations in the least-loaded class can increase over time owing to Muller’s ratchet. To keep constant the average number of mutations in the eDNA pool relative to the number of mutations in the least-loaded class (so that the dynamics is stationary over time), the fitness in was normalized by the fitness of the least leaded class to obtain the expression for εi shown previously. Note that assuming εi is more unfavorable to the hypothesized advantage of HGT than assuming because is nearly independent of i when .
In addition, the model assumed spontaneous decay of eDNA at rate per cell generation. The eDNA decay rate determines the turnover of the eDNA pool, i.e., how synchronous the genetic content of the eDNA pool is relative to that of the population (e.g., when , the entire eDNA pool is replaced with the new input every generation). Although the absolute concentration of eDNA is likely to influence the frequency of HGT (Levin and Cornejo 2009), in the model this influence was assumed to be absorbed by r so that the value of deDNA does not influence the frequency of HGT. To indicate this fact, deDNA will be henceforth referred to as the rate of eDNA turnover.
With the aforementioned assumptions, the dynamics of the eDNA pool was simulated as follows. The state of the eDNA pool was described by the vector where eij is the copy number of allele j (i.e., j mutations) from locus i present in the eDNA pool. In each generation, reduced by deDNAeij (if deDNAeij > 1, it was normalized to the nearest integer; otherwise, the decrease of eij was drawn from a binomial distribution with the number of trials eij and success probability deDNA). Then, N cells were randomly chosen with replacement from the population of the current generation with the probability εi. All alleles in the genomes of the chosen cells were added to the eDNA pool, leading to the increase of eij.
The simplest version of the model assumed that the population is genetically isolated from other populations. However, to investigate the effect of genetic exchange between semi-isolated populations, we also considered an extended model in which the population was subdivided into 16 subpopulations of an equal size. The subpopulations were arranged in a two-dimensional square grid (4 × 4 subpopulations) with toroidal boundaries. Each subpopulation exchanged individuals with the four nearest-neighbor subpopulations via random migration. The number of migrants between each pair of neighboring subpopulations was drawn from the binomial distribution with the success probability Dpop and Ns/4 trials where Ns is the size of a subpopulation (Ns = N/16). Each subpopulation had its own pool of eDNA that exchanged alleles with the four nearest-neighbor pools via diffusion. The copy number of each allele diffusing out of a pool was defined as a fraction DeDNA/4 per direction per generation (if DeDNAeij/4 > 1, it was normalized to the nearest integer; otherwise, the numbers of copies diffusing out were drawn from a multinomial distribution with the number of trials eij and success probability DeDNA/4 per direction).
Organization of the simulation
The dynamics of the model in one generation was simulated in four steps, in the following order:
1. The Wright-Fisher process. In this step, the population of the current generation died and was replaced by that of the next generation.
2. Production and decay of the eDNA pools.
3. Diffusion of both the microbial populations and eDNA.
4. Mutation and HGT.
Characterization of Muller’s ratchet
With HGT disabled, the model is known to have the following family of stationary distributions that are distinguished by the number of deleterious mutations in the least-loaded class : , ⋯, , where k is a non-negative integer (Haigh 1978). A finite-size population cannot stay in any of these distributions indefinitely because of stochasticity. The absence of back mutations implies that mLLC increases over time, hence the accumulation of deleterious mutations (Haigh 1978). In general, such accumulation can result from two different processes (Charlesworth et al. 1993): the fixation of deleterious mutations via genetic drift or the operation of Muller’s ratchet per se, i.e., the extinction of the least-loaded class as the result of stochastic fluctuations. Muller’s ratchet can be reversed by recombination only if a mutation has not been fixed (unless eDNA turnover is very slow, as described later in this article). Thus, it is beneficial to know through which of these processes the accumulation of mutations occurs in the model. Thus, we first analyzed the model in the absence of HGT (i.e., r = 0). The results of the analysis essentially confirmed the conclusion of Charlesworth and Charlesworth (1997) with additional data. Although these results in part reproduce the work of Charlesworth and Charlesworth (1997), we present them here because they serve as the baseline with which to compare the results presented in the later sections.
The simulations show that the tempo and mode of mutation accumulation in the model depend on whether the value of () is much greater than unity or not, as expected from previous work (Rouzine et al. 2008; Neher and Shraiman 2012). For (slow accumulation regime), mLLC and the number of mutations fixed in the population (mfix) increase over time in a very similar fashion (Figure 1A), so that the number of segregating mutations in the least-loaded class (mLLC − mfix) is at most one. Each increase in mLLC is followed by an increase in mfix with a delay of approximately a few thousand generations before the next increase occurs to mLLC (Figure 1B). Thus, the extinction of the least-loaded class (i.e., a click of Muller’s ratchet) always precedes the fixation of a mutation (Charlesworth and Charlesworth 1997).
The dynamics of mutation accumulation was further analyzed by measuring the average gene diversity (“heterozygosity”) of the least-loaded class (denoted ) and that of one but the least-loaded class (; Nei 1987). is defined as where is the frequency of allele j in locus i in the least-loaded class (Nei 1987). is defined in a similar manner. The results show that and display the following dynamics. When Muller’s ratchet clicks (i.e., when the value of mLLC increases), one but the least-loaded class becomes the new least-loaded class by definition, so that the value of jumps to the value of immediately prior to the click of the ratchet (Figure 1C). This value of is always nearly 0.02, indicating that, immediately prior to a click of Muller’s ratchet, one but the least-loaded class consists of l genotypes (one mutation in one of the l loci), whose frequencies roughly equal to each other ( for ). Thus, immediately after the ratchet clicks, the new least-loaded class consists of l genotypes. Because these genotypes have an equal fitness value, the gene diversity decreases to zero over time owing to genetic drift (Figure 1C). Almost simultaneously with hitting zero, reaches a steady state value of approximately 0.02 mentioned previously (Figure 1C), and the value of mfix increases, indicating the fixation of a deleterious mutation (Figure 1B). In summary, the accumulation of mutations occurs as a cycle of three distinct steps (Figure 2):
1. Extinction of the least-loaded class due to stochasticity (i.e., Muller’s ratchet);
2. Gradual decline of the genetic diversity in the new least-loaded class due to genetic drift;
3. Eventual fixation of a deleterious mutation in the entire population.
For , mutations accumulate more rapidly than when (Figure 1D). The value of mLLC is incremented multiple times before the value of mfix reaches the value of mLLC (Figure 1E). This result indicates that, unlike the case of , cycles of mutation accumulation can overlap with each other in time: a new cycle can begin with a click of Muller’s ratchet (i.e., step 1 in the previous paragraph) before the previous cycle undergoes the fixation of a mutation (i.e., step 3). Moreover, the dynamics of and are no more synchronized, in that the decrease of after a click of Muller’s ratchet is delayed with respect to the decrease of (cf. Neher and Shraiman 2012). Likewise, the increase of mfix also is delayed with respect to the decrease of . These results indicate that the decline of gene diversity after a click of Muller’s ratchet (i.e., step 2) is a sequential process, in that gene diversity is consecutively lost in different genotype classes distinguished by the number of mutations (i.e., first, decreases to zero, and then decreases to a steady-state level, and so on). Taken together, the results indicate that the three-step picture of mutation accumulation described in the previous paragraph remains essentially valid in the rapid ratchet regime () as well.
Therefore, Muller’s ratchet is the dominant process leading to the accumulation of mutations in the model with HGT disabled (rather than the fixation by drift of mutations; Charlesworth and Charlesworth 1997). However, a click of Muller’s ratchet eventually results in the fixation of a mutation, so that the number of fixed mutations (mfix) far exceeds the number of segregating mutations in the least-loaded class (mLLC − mfix). This situation is in stark contrast with the result with a diploid, sexual model that assumes no recombination and weak dominance, which shows that the number of fixed mutations is much smaller than the number of segregating mutations (Charlesworth et al. 1993). The most probable cause of this discrepancy is the prevention of mutation fixation in diploid populations due to homozygote disadvantage, which is absent in our haploid model. The resulting paucity of segregating mutations in the haploid model seems to work against the putative advantage of HGT; later we show how this limitation might be overcome by population subdivision.
Effect of HGT on Muller’s ratchet
We next examined the effect of HGT on Muller’s ratchet. For simplicity, a rapid turnover of eDNA was assumed by setting deDNA = 1 (the effect of slow eDNA turnover is presented in the next section). In this case, assuming a very large population, the average number of mutations in the eDNA pool per l loci (i.e., per genome) can be approximated by U/s + 1, where U/s is the average number of mutations in the population at the steady state with mLLC = 0 (Redfield et al. 1997; File S1). Because the mean number of mutations in the eDNA pool is greater than that in the population, HGT on average introduces more deleterious mutations into the population than it removes and so decreases the steady-state population size of the least-loaded class.
To assess the net effect of HGT on Muller’s ratchet, we measured the rate of mutation accumulation Δmfix /Δt as a function of r (HGT rate) for different combinations of N and values. The results show that HGT prevents the accumulation of mutations (i.e., Δmfix /Δt = 0) if its frequency is sufficiently high (Figure 3). Moreover, the frequency of HGT required to stop Muller’s ratchet decreases with the increase of N, indicating that the preventive effect of HGT increases with N. This trend is seen even in the slow ratchet regime (), in which mutation accumulation accelerates with the increase of N in the absence of HGT (Figure 3) (note that because and s are fixed, U was set to a greater value for a greater value of N; thus, mutation accumulation can accelerate with the increase of N).
To understand why HGT can prevent Muller’s ratchet despite its average deleterious effect due to the high mutation load of eDNA that comes from dead cells, we considered the following simple Wright-Fisher model. The model assumed two populations: the population of the least-loaded class n0 and that of the other classes n1 (the model is thus equivalent to one-locus, two-allele, haploid model). The total population size N = n0 + n1 was constant. The value of n0 at time t + 1 was drawn from a binomial distribution with N trials and success probability p′ defined as follows:where p = n0/N and () at time t, s′ is the fitness effect of mutation, U is the mutation rate from n0 to n1 (back mutations ignored), is the average fitness, and r01 and r10 are the rates of HGT. For simplicity, HGT was assumed to convert n0 into n1 at rate r01 and, conversely, n1 into n0 at rate r10. The effect of eDNA was implicitly incorporated into the model by the choice of parameters r01 and r10 as follows. If r01 and r10 satisfy the following condition, HGT decreases the steady-state level of n0 compared with the steady-state level reached in the absence of HGT (i.e., r01 = r10 = 0):(the left-hand side of this inequality is the per-generation conversion of n1 into n0, whereas the right-hand side is that of n0 into n1, when n0 and n1 are set to the steady-state level achieved in the absence of HGT).
With this simplified model, we asked whether HGT facilitated or antagonized the survival of n0. The mean time to extinction of n0, which is a proxy for Δmfix/Δt in the full model, was measured as a function of r01 and r10 for different combinations of N and U′ values. The results show that HGT increases the mean time to extinction of n0 whether or not HGT decreases the steady state level of n0, as long as the per-generation production of the least-loaded class by HGT () is on the order of unity or greater (Figure 4; note that for the parameters used in the simulations). This result makes intuitive sense: the least-loaded class cannot go extinct if it is continually produced by HGT in every generation; and as long as such production is guaranteed, the steady-state population size of n0 is immaterial to the extinction of n0. To check whether this explanation also applies to the full model, we investigated the consistency between the two models from a few different angles as described below.
First, the simplified model indicates that the required frequency of HGT to prevent Muller’s ratchet decreases as N increases because increasing N increases the production rate of the least-loaded class (with kept constant). The same result is produced by the full model as already described previously.
Next, Figure 4 shows that whether HGT increases the time to extinction depends on the value of r10 much more strongly than on the value of r01 (as seen from the fact that the boundary between yellow and red regions is nearly horizontal). We thus hypothesized that, in the full model, HGT prevents Muller’s ratchet most effectively when the probability of HGT producing the least-loaded class is maximized. This probability can be approximated by the probability P10 that an HGT event transforms one but the least-loaded class () into the least-loaded class (), because the contributions from other mutant classes likely involve more than one HGT event. P10 can be calculated as under the assumption that N is very large (File S1). P10 takes the maximum value when ( when ; File S1). The non-monotonic dependency of P10 on l can be intuitively understood as follows: decreasing l can decrease P10 because it increases the probability that an allele randomly drawn from the eDNA pool contains a mutation that is not carried by the least-loaded class (if l = 1, this probability is one). However, increasing l also can decrease P10 because it decreases the probability that HGT occurs to the very locus in which one but the least-loaded class has the mutation not carried by the least-loaded class. Therefore, P10 takes the maximum value at an intermediate value of l. By contrast, the probability P01 that HGT introduces one or more mutations in the least-loaded class is (File S1), which is a monotonically decreasing function of l because increasing l decreases the probability that an allele randomly drawn from the eDNA pool contains mutations. Thus, the aforementioned hypothesis can be reformulated as follows: Δmfix/Δt reaches the minimum value when l assumes the value . To examine this hypothesis, we measured Δmfix/Δt as a function of l for various combinations of N and (Figure 5). The results show that Δmfix/Δt takes the minimum value when the value of l is on the same order of magnitude as (around 10 rather than 100), in accord with the hypothesis.
Finally, the simplified model shows that is sufficient for HGT to prevent the extinction of . To check if a similar condition applies to the full model, we calculated the per-generation production of the least-loaded class by HGT occurring to one but the least-loaded class (). For , sharply drops at for N = 106 and 105 (Figure 3). The values of in these cases are approximately 0.6 and 0.4 for N = 106 and 105, respectively, in qualitative agreement with the expectation that be on the order of unity (note that is expected because the contributions from the more-loaded classes are ignored and is a stronger condition than ). For , sharply drops at for N = 106 and 105. The values of in these cases are approximately 0.8 and 0.6 for N = 106 or 105, respectively, again in qualitative agreement with the expectation. (The results were also consistent for smaller values of N, at least in the orders of magnitude.)
Taken together, the consistency of the results indicates that the simplified model captures, at least qualitatively, the mechanism by which HGT prevents Muller’s ratchet in the full model. Namely, HGT can prevent Muller’s ratchet even if it decreases the steady-state population size of the least-loaded class because it can hinder the extinction of the least-loaded class by continually generating it from the more-loaded classes.
Effect of the eDNA turnover rate on the prevention of Muller’s ratchet by HGT
We further investigated the effect of the eDNA turnover rate on the ability of HGT to prevent Muller’s ratchet. As described previously, the accumulation of mutations in the model is initiated by a click of Muller’s ratchet followed by the fixation of a mutation. Once a mutation is fixed (i.e., once a good allele is lost), HGT cannot undo it because the only source of HGT in the model is eDNA derived from the members of the same population. However, slow turnover of eDNA can delay the loss of good alleles from the eDNA pool and thereby might help HGT prevent the accumulation of mutations.
To examine this possibility, we measured as a function of deDNA for various combinations of N and values (Figure 6). The results show that Δmfix/Δt decreases with the decrease of deDNA. However, for this effect to be significant, eDNA turnover has to be slower than the population turnover (i.e., generation time) by more than three orders of magnitude (deDNA < 10−3), a condition that appears to be unrealistic in natural environments (Nielsen et al. 2007; Pietramellara et al. 2009).
To understand why such a severe condition applied, we measured the “quality” of eDNA defined as where xki is the fraction of the least-loaded class having allele i (i.e., i mutations) in locus k, and is the fraction of allele j from locus k in the eDNA pool. qeDNA indicates the potential of the eDNA pool to remove deleterious mutations in the genomes of the least-loaded class (e.g., qeDNA is zero if the eDNA pool contains no allele that can remove mutations in the least-loaded class). qeDNA was measured in the absence of HGT for various values of deDNA. The results show that, in the slow ratchet regime (), the dynamics of qeDNA closely follows that of when eDNA turnover is not very slow (deDNA ≥ 10−2) (Figure 1C and Figure 7, A and B). When eDNA turnover is very slow (deDNA ≤ 10−3), the decay of qeDNA lags behind the decline of (Figure 7C). This result indicates that eDNA turnover has to be slower than the decline of to delay the decay of qeDNA significantly. Because the decrease of is attributed to genetic drift as mentioned previously, its timescale is determined by the population size of the least-loaded class, which was approximately 1000 in the simulations for the slow ratchet regime ( and ). Thus, deDNA must be <10−3 to cause any significant effect.
The aforementioned argument might imply that the effect of slow eDNA turnover could be more significant if the population size of the least-loaded class is small. However, this turns out not to be the case. In the fast ratchet regime ( and ), the dynamics of qeDNA no more closely follows the dynamics of , and it still takes more than 500 generations for qeDNA to decrease to zero when eDNA turnover is rapid (Figure 7, D−F). Therefore, deDNA still has to be smaller than 1/500 to have any significant effect (File S1 contains a more detailed explanation of the results described in this and the previous paragraphs).
In summary, the fixation of mutations is an evolutionary process that is driven by selection and drift and thus occurs on relatively long timescales compared with the generation time of individuals, so that eDNA turnover has to be commensurately slow to have any appreciable effect.
Effect of population subdivision on the prevention of Muller’s ratchet by HGT
Next, we considered population subdivision. On the one hand, population subdivision can accelerate accumulation of mutations because it enhances the effect of genetic drift as shown in previous studies (Higgins and Lynch 2001; Combadao et al. 2007). On the other hand, subdivision might help HGT prevent Muller’s ratchet as follows. Let us suppose that the isolation between subpopulations is so strong that the accumulation of mutations in different subpopulations is independent of each other in the absence of HGT. In this case, a mutation must be fixed in every subpopulation independently before being fixed in the entire population. This situation is likely to lead to an increased number of segregating mutations when the entire population is considered (rather than when each subpopulation is considered separately). The increased number of segregating mutations by itself does not affect accumulation of mutations if recombination requires direct contact between individuals as in many animals. HGT, however, does not require such direct contact. If subpopulations share a common eDNA pool because of the rapid transport of DNA molecules between their habitats, increasing the number of segregating mutations is expected to increase the number of mutations that are introduced or removed by HGT events. This should enhance the preventive effect of HGT on Muller’s ratchet because removal of mutations by HGT is more important than introduction of mutations by HGT according to the simplified Wright-Fisher model presented above. The key question is whether this potential advantage of population subdivision can outweigh the disadvantage of population subdivision associated with the enhanced genetic drift.
To address the aforementioned question, we measured as a function of population migration rate (Dpop) in a model incorporating population subdivision (Materials and Methods). was defined as the average of mfix measured within each subpopulation. For simplicity, all subpopulations were assumed to share a common eDNA pool (i.e., DeDNA = ∞; we will consider finite DeDNA later). The results show that, in the absence of HGT (r = 0), the system becomes more prone to Muller’s ratchet (i.e., increases) as Dpop decreases (Figure 8, A and B). This outcome is expected, given the disadvantage of population subdivision due to enhanced genetic drift. By contrast, in the presence of HGT (r ≥ 10−4), the dependency of on Dpop is non-monotonic: the system initially becomes more prone to Muller’s ratchet as Dpop decreases, but becomes less so as Dpop decreases further, with the position of the maximum along the Dpop axis shifting as a function of r (Figure 8, A and B). The comparison against the curve for r = 0 shows that the reduction in caused by HGT increases as Dpop decreases, indicating that population subdivision increases the effect of HGT. At a high HGT rate (r = 0.01), in the fast ratchet regime (), the system is more resistant to Muller’s ratchet at lower values of Dpop (≤10−5) than at greater values of Dpop (≥10−1) for which the system can be considered well-mixed (Figure 8B). In the slow ratchet regime (), however, this situation cannot be observed because falls to zero and so shows no difference between low and high values of Dpop within the accuracy of the simulations (Figure 8A). In addition to the aforementioned results, was measured as a function of Dpop for other combinations of N and values; the results were qualitatively the same as described previously (File S1). Taken together, these results indicate that the advantage of population subdivision can outweigh the associated disadvantage, provided that the frequency of HGT is sufficiently high, and subpopulations share a common eDNA pool.
To elucidate the mechanism by which population subdivision helps HGT prevent Muller’s ratchet, we first sought to test the expectation that population subdivision increases the number of segregating mutations at the level of the entire population. To this end, we measured the average gene diversity of the least-loaded class within and between subpopulations (Nei 1987). The average gene diversity within one subpopulation is defined as where xyki is the frequency of allele i in locus k in the least-loaded class of subpopulation y (the least-loaded class was defined separately for each subpopulation). The average gene diversity within subpopulations HS is defined as the average HSy over all subpopulations (weighted by the population size of the least-loaded class in each subpopulation). The average gene diversity in the entire population HT is defined similarly with xyki replaced by the frequency of an allele in the least-loaded classes from all subpopulations combined together. The average gene diversity between subpopulations was defined as HT − HS (Nei 1987). HT and HS were measured in the absence of HGT. Simulations show that the values of HT and HS initially increase and then saturate over time (results not shown). The time averages of HT and HS after or near saturation were obtained as a function of Dpop (Figure 9A). The result shows that HT increases as Dpop decreases, whereas HS remains close to zero for the entire range of Dpop. This result indicates that population subdivision increases the number of loci that are polymorphic between subpopulations, whereas within subpopulations almost all loci remain monomorphic. Therefore, population subdivision increases the number of segregating mutations at the level of the entire population as expected.
The increase in the number of segregating mutations is expected to drive an increase in the number of mutations that are introduced or removed by HGT events. To test this expectation, we measured the quality of the eDNA pool defined above (qeDNA) and the “toxicity” of the eDNA pool defined as (qeDNA and teDNA were calculated for each subpopulation and averaged over subpopulations using the population sizes of the least-loaded classes as weights). teDNA indicates the potential of the eDNA pool to introduce deleterious mutations in the genomes of the least-loaded class. Note that the difference teDNA − qeDNA is equal to the average number of mutations introduced (or removed if the value is negative) by one HGT event occurring to an individual of the least-loaded classes. Thus, − qeDNA can be interpreted as the contribution to the average number of mutations introduced by HGT by those HGT events that decrease the number of mutations, whereas teDNA as the contribution by those events that increase the number of mutations. Simulations show that the values of qeDNA and teDNA increase and then saturate over time except for Dpop = 0 (result not shown; for Dpop = 0, saturation is likely to be reached, but probably requires an exceedingly long time). The time averages of qeDNA and teDNA after saturation were obtained as a function of Dpop in the absence of HGT (Figure 9B; for Dpop = 0, the values reached toward the end of the simulation are shown, thus indicating a lower bound). The result shows that population subdivision increases the number of mutations removed by an HGT event (qeDNA) as well as the number of mutations introduced by an HGT event (teDNA), thus confirming the expectation (Figure 9B). Moreover, it leaves nearly constant the average (i.e., net) number of mutations introduced by an HGT event (teDNA − qeDNA). In other words, population subdivision increases the conversion of the least-loaded class by HGT in the direction of both increasing and decreasing the number of mutations without changing the magnitude and direction of the net conversion. According to the simplified Wright-Fisher model described before, removal of mutations by HGT (r10) has a greater influence on the extinction of the least-loaded class than introduction of mutations by HGT (r01). Although r10 and r01 are concerned with HGT occurring in different genotypes, they are analogous to qeDNA and teDNA, respectively. This analogy points to the advantage of increasing qeDNA and teDNA with teDNA − qeDNA kept constant to prevent Muller’s ratchet.
Whether population subdivision is advantageous or disadvantageous depends on the combination of parameters. To examine what appears to be the most important aspect of this parameter dependence, we next considered the finite diffusion of eDNA. To this end, the value of was compared between Dpop = 10−5 (subdivided population) and Dpop = 10−1 (well-mixed population) for various combinations of deDNA and DeDNA values. As expected, the results show that is smaller for Dpop = 10−5 than for Dpop = 10−1 when deDNA is sufficiently small or DeDNA is sufficiently large (Figure 10). The values of DeDNA and deDNA demarcating the region of the parameter space for which population subdivision is advantageous do not appear unrealistic, although it is difficult to judge whether the conditions represented by this region are achievable in natural environments (see also Discussion).
The results presented here suggest that HGT can prevent the operation of Muller’s ratchet in prokaryotic populations, even if on average HGT introduces more deleterious mutations than it removes. The avoidance of Muller’s ratchet via transformation and recombination with eDNA might appear somewhat paradoxical because on average eDNA is expected to carry more deleterious mutations than the DNA in live prokaryotic cells. Indeed, it has been long recognized that “sex with dead cells” is a dubious proposition (Redfield 1988, 2001; Redfield et al. 1997). However, the modeling results indicate that, even though transformation via eDNA is expected to increase the mean mutation load and hence decrease the average fitness of a population, it nevertheless can stop Muller’s ratchet. This appears to be the case because HGT provides for the chance to eliminate deleterious mutations, leading to the continual restoration of the least-loaded class. In other words, HGT prolongs the persistence of the least-loaded class in the population in the face of stochastic fluctuations, hence the prevention of Muller’s ratchet.
The efficacy of HGT in the prevention of Muller’s ratchet depends on the stability of eDNA, its diffusion rate, and population subdivision. The contribution of eDNA stability (i.e., slow turnover) is intuitively clear because it delays elimination of high-fitness alleles from the population and so increases the chance that these alleles offset Muller’s ratchet via HGT. However, the simulations show that, for this effect to be substantial, the characteristic lifetime of eDNA should exceed the generation time by orders of magnitude (Figure 6), suggesting that eDNA stability alone is not a key factor in the evolution of microbial populations. If, however, a population is subdivided, enhanced eDNA stability facilitates the transport of eDNA across subpopulations and thus contributes to the ability of HGT to prevent Muller’s ratchet. Under these conditions, for the effect to be significant, the characteristic lifetime of eDNA does not have to be exceedingly long, depending on the diffusion rate of eDNA (Figure 10).
If the transport of eDNA is sufficiently fast, HGT can prevent Muller’s ratchet more efficiently in a subdivided population than in an undivided population. As a result, a subdivided population can be more resistant to Muller’s ratchet than an undivided population of an equal overall size despite the fact that population subdivision accelerates Muller’s ratchet in the absence of HGT owing to the enhanced effect of genetic drift. Put differently, to maintain genomic information in the face of deleterious mutations and genetic drift, it is advantageous to partition individuals into multiple (smaller) subpopulations and let them “cross-reference” each other’s genetic information through recombination rather than collect all individuals in one population and thereby maximize the efficacy of natural selection. Note, however, that although population subdivision can be beneficial, it is not necessary for HGT to prevent Muller’s ratchet.
How do these findings relate to the previous population genetics studies on the advantage of recombination? Previous studies have shown that recombination is beneficial (compared with no recombination) when a population suffers from the reduced efficacy of natural selection because of negative linkage disequilibrium (i.e., biased associations on chromosomes between beneficial and deleterious alleles) which arises from genetic drift (the Hill-Robertson effect) or synergetic epistasis (Hill and Robertson 1966; Felsenstein 1974; Feldman et al. 1980; Kondrashov 1988, 1993; Otto and Lenormand 2002). Our present results are consistent with these findings because they show that the reduction in the rate of mutation accumulation caused by HGT is greater in a subdivided population than in an undivided population (Figure 8, A and B). Indeed, population subdivision causes negative linkage disequilibrium in many pairs of loci in the least-loaded classes at the entire population level because the least-loaded classes in different subpopulations have mutations in different loci (this can be seen from the fact that HT > 0 and HS ≈ 0 when Dpop ≈ 0). In addition, the results show that the rate of mutation accumulation itself can be lower in a subdivided population than in an undivided population, despite the disadvantage of population subdivision due to enhanced drift. This conclusion does not immediately follow from the previous findings because the fact that negative linkage disequilibrium increases the benefit of recombination (compared with no recombination) does not necessarily mean that negative linkage disequilibrium is beneficial in preventing the accumulation of mutations in the presence of recombination.
Evolutionary maintenance of genomic information based on population subdivision and recombination hinges on the fact that genetic exchange between cells happens via a shared pool of eDNA. On the one hand, this implies that the indirectness of HGT might actually be an asset rather than a drawback for prokaryotes because HGT provides the possibility of exchange between physically separated subpopulations that accumulate different mutations.
On the other hand, indirect recombination imposes strict conditions under which population subdivision is advantageous in terms of Muller’s ratchet, requiring restricted migration between populations and rapid transport of eDNA molecules. Although at present it is difficult to judge whether such conditions are fulfilled in natural environments, the available data seem to suggest soils as an environment that is conducive to these processes. The distribution of bacteria in natural soils is “patchy,” i.e., bacteria typically occur as small colonies separated from each other, and their locations are restricted by the availability of soil micro-pores, water, and organic substances (Foster 1988; Stotzky 1989; Grundmann 2004; Young et al. 2008; Vos et al. 2013). Moreover, bacterial motility can be highly restricted by water availability (Dechesne et al. 2010). In addition, soil bacteria display high levels of genetic diversity within local populations (Grundmann 2004). For example, Agrobacterium biovar 1 and Nitrobacter-like bacteria within minute soil samples (<1 mm in diameter) display genetic divergence as great as that between strains sampled from different geographical areas (Grundmann and Normand 2000; Vogel et al. 2003, although Agrobacterium biovar 1 is not known to be naturally competent, this bacterium readily undergoes homologous recombination; Costechareyre et al. 2009). Also, genomes of Pseudomonas stutzeri display extremely high gene diversity (HT > 0.8) and high frequency of null alleles (>80% of the examined strains fail to exhibit the activity of at least one of the 20 enzymes examined; Rius et al. 2001). Taken together, these findings suggest that bacterial populations in unsaturated soils are strongly subdivided and genetically heterogeneous within local environments. This population structure is compatible with a situation in which HGT and population subdivision jointly prevent the operation of Muller’s ratchet (if eDNA transport is sufficiently rapid as well).
Moreover, eDNA in soils becomes highly stable against hydrolysis by nucleases when it is absorbed to clay minerals, sand particles, and humic acids without losing the ability to transform competent bacteria (Vries and Wackernagel 2005; Nielsen et al. 2007; Pietramellara et al. 2009). Natural transformation assays indicate that DNA added to nonsterile soils retains transformation ability for 3−15 days (Gallori et al. 1994; Sikorski et al. 1998). In addition, eDNA can be transported within unsaturated soils through water capillarity or leaching (Ceccherini et al. 2007; Pote et al. 2010). Although the difference between bacteria and eDNA in terms of transport efficiency within soils remains to be evaluated, the obvious size difference suggests that eDNA is more readily transported than bacteria. Taken together, these findings suggest soil prokaryote systems as a prime candidate for a system in which HGT and population subdivision could be important for the long-term maintenance of genomic information in the face of Muller’s ratchet. Furthermore, it appears plausible that both the molecular machinery of natural transformation and GTAs evolved, at least in part, as adaptations facilitating HGT and overcoming the limitations on HGT imposed by the aforementioned conditions, e.g., GTAs can increase the stability of DNA emitted from the cell.
Is Muller’s ratchet an important factor of evolution in prokaryotes? Its role can be questioned because many bacteria have large population sizes in which the extent of stochasticity is limited and that are accordingly immune to Muller’s ratchet (but see also the discussion of soil bacteria above in this section). However, bacteria are also known to go through frequent population bottlenecks under stress conditions (Levin et al. 2000; Aertsen and Michiels 2004; Wolf et al. 2005), and it is under these conditions that Muller’s ratchet and mechanisms that can overcome it could become relevant. Indeed, competence has been shown to increase stress resistance in bacteria (Engelmoer and Rozen 2011) and induction of the competence regulons is an important part of bacterial stress response (Claverys et al. 2006; Prudhomme et al. 2006). Although the biology of the GTAs is not nearly as well understood as that of competence and transformation, stress-induced prophage replication that is likely to potentiate HGT has been demonstrated (Garcia-Russell et al. 2009).
To conclude, the results of the modeling study described here indicate that HGT can prevent the operation of Muller’s ratchet, and this effect can be enhanced by population subdivision. Thus, through preventing the action of Muller’s ratchet, HGT could be a necessary condition for the long-term survival of prokaryotic populations. The benefits of HGT for prokaryotes certainly are not limited to the prevention of Muller’s ratchet and additionally include acquisition of novel genes and functions (see, e.g., Levin and Cornejo 2009; Wylie et al. 2010). Jointly, these advantages could be the driving forces behind the evolution of dedicated mechanisms for HGT.
We thank Paulien Hogeweg, Alexander Lobkovsky, Nen Saito, and Yuri Wolf for helpful discussions. This research was in part supported by the Intramural Research Program of the National Institutes of Health, National Library of Medicine, by the Japan Society for the Promotion of Science Research Fellowship for Japanese Biomedical and Behavioral Researchers at the National Institutes of Health, and by the Japan Society for the Promotion of Science Research Fellowship for Young Scientists.
Supporting information is available online at http://www.g3journal.org/lookup/suppl/doi:10.1534/g3.113.009845/-/DC1
Communicating editor: B. J. Andrews
- Received September 4, 2013.
- Accepted December 12, 2013.
- Copyright © 2014 Takeuchi et al.
This is an open-access article distributed under the terms of the Creative Commons Attribution Unported License (http://creativecommons.org/licenses/by/3.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.