Selection against Heteroplasmy Explains the Evolution of Uniparental Inheritance of Mitochondria

Summary

Mitochondria contain genes that encode the machinery needed to power cells. Unlike the nuclear genome, the mitochondrial genome is typically inherited from one parent only (uniparental inheritance). The most common explanation for uniparental inheritance is the genomic conflict theory, which states that uniparental inheritance evolved to prevent the spread of ‘selfish’ mitochondria that replicate quickly but produce energy inefficiently. Current explanations have a major problem: when using realistic parameters, mathematical models cannot show that uniparental inheritance can replace biparental inheritance. Clearly, we need a new explanation that fits with standard population-genetic theory. Recent evidence suggests cells may incur a cost when they carry multiple types of mitochondria. Here we show mathematically that uniparental inheritance could have evolved to avoid the costs of maintaining multiple mitochondrial lineages within a cell. Our results explain the long-standing evolutionary mystery of uniparental inheritance and provide insight into the evolution of mating types and binary sexes. Selection against heteroplasmy also has implications for the evolution of the mitochondrial genome because new mitochondrial haplotypes always lead to heteroplasmy before becoming fixed in the population. Thus, selection against heteroplasmy may explain why mtDNA coding-genes have slower substitution rates than analogous genes within the nucleus.

Introduction

During sexual reproduction, offspring receive two genomes: nuclear genomes from both parents and haploid cytoplasmic genomes, contained in mitochondria and chloroplasts (in plants and algae), usually from one parent. Although uniparental inheritance is nearly ubiquitous, the reasons behind its evolution remain unresolved [1, 2]. Cells contain multiple mitochondria, and the mitochondrial genome (mtDNA) encodes polypeptide subunits of the electron transport chain, which the cell uses to generate ATP via oxidative phosphorylation [2]. If mutations increase mtDNA replication rate but simultaneously decrease respiration, then increased mtDNA fitness comes at the expense of cell and organism fitness [3–5]. Nuclear and mitochondrial genomes are thus potentially in conflict. The genomic (or selfish) conflict theory argues that uniparental inheritance evolved because biparental inheritance facilitates the spread of such selfish mitochondria [1, 3–6]. Although the conflict theory has been the predominant explanation for uniparental inheritance for over three decades [3, 4], other explanations exist. A second theory suggests that uniparental inheritance facilitates the removal of deleterious mutations. Uniparental inheritance decreases variation of mtDNA within cells, but increases variation between cells, allowing purifying selection against cells with increased mutation load [1, 7]. A third hypothesis argues that because the oxidative phosphorylation pathway is composed of interacting nuclear- and mitochondrial-encoded polypeptides, uniparental inheritance optimizes mitochondrial-nuclear coadaptation by maintaining coevolved mitochondrial-nuclear combinations [1, 8]. While uniparental inheritance spreads in mathematical models of the above hypotheses [1, 5, 6], it cannot replace biparental inheritance under realistic assumptions and parameter values [1, 5]. Thus, despite decades of theoretical work, we still lack a convincing explanation for why uniparental inheritance is widespread amongst extant organisms [1, 2].

Although uniparental inheritance is the general rule in eukaryotes, there are a few exceptions. Probably the best-known exception is baker’s yeast (Saccharomyces cerevisiae) in which both parents contribute mitochondria to offspring [9, 10]. However, the repeated division of cells that contain two mitochondrial lineages (heteroplasmy) leads to cells that contain a single type of mitochondria (homoplasmy) [9, 10]. Another example is the male bivalve (Mytilus), which also inherits mitochondria from both parents. But in this case maternal and paternal mitochondria do not mix within single cells, as maternal mitochondria segregate to the soma while paternal mitochondria segregate to the gonads [11]. Thus, even when mitochondria are inherited from both parents, heteroplasmy is avoided. Recent experimental evidence suggests that this is because heteroplasmy imposes a cost on the organism. A study on mice found that the mere mixing of different, but phenotypically normal, mitochondria within a cell leads to physiological and behavioral abnormalities [12]. Could uniparental inheritance have evolved simply because carrying multiple mitochondrial types imposes a cost on the organism? Here we use a mathematical model to explore whether selection against heteroplasmy could have led to the evolution of uniparental inheritance.

Basic model description

Our model is based on an idealized life cycle of a single-cell diploid eukaryotic organism, such as the algae Chlamydomonas reinhardtii. Diploid cells contain n mitochondria and haploid cells have n/2 mitochondria. All mitochondria are initially wild type but mitochondria can mutate from wild type to mutant (and vice versa). The starting population contains haploid gametes with a nuclear allele regulating biparental inheritance (B). Gametes are evenly split between two nuclear self-incompatible mating types (B₁ and B₂). In the basic model, we assume no recombination between the mitochondrial inheritance and mating type loci because these are tightly linked in many isogamous organisms [9] (later we explore recombination and no mating types). Cell types are characterized by the proportion of wild type and mutant mitochondria that they carry and their nuclear allele (haploid) or genotype (diploid).

Our life cycle has four discrete stages and is similar to the life cycles used in previous models [1, 5, 8]. Since we begin with a population of gametes, the first stage is random mating. Here, gametes randomly mate with the opposite mating type to produce diploid cells. Matings are controlled by the nuclear allele in gametes. In biparental inheritance (between B₁ and B₂ gametes), both gametes contribute mitochondria to the B₁B₂ diploid cells (see later for uniparental inheritance). The second stage is mutation. Each mitochondrion can mutate to the other haplotype with probability μ. The third stage is selection. Here, diploid cells have a relative fitness based on the proportion of each haplotype in the cell. We assume that fitness decreases as the level of heteroplasmy increases. The fourth stage is meiosis, where diploid cells produce gametes that contain a single nuclear allele and n/2 mitochondria. As mitochondria are stochastically partitioned into gametes [9], diploid heteroplasmic cells produce gametes with varying degrees of heteroplasmy.

First, we let the population of B₁ and B₂ gametes reach mutation-selection equilibrium. We then simulate a mutation leading to uniparental inheritance of mitochondria by converting a small proportion (10⁻²) of B₁ gametes to U₁ gametes. We assume no further mutations between B and U alleles. Matings between U₁ and B₂ gametes result in uniparental inheritance, in which the U₁B₂ cell inherits mitochondria from U₁ alone. (Matings between U₁ and B₁ are not possible as they are the same mating type.) The population now consists of three alleles (U₁, B₁ and B₂) and two genotypes (U₁B₂ and B₁B₂). The model tracks the proportion of each cell type at each stage of the life cycle. U₁ spreads at the expense of B₁ when uniparental inheritance is more advantageous than biparental inheritance (the frequency of B₂ always remains at 0.5), and the simulation ends when the alleles reach equilibrium (see Model and S1–S6 Model for details of the model).

To explore whether a cost to heteroplasmy could have led to the evolution of uniparental inheritance, we study several scenarios. We first examine the simplest case, where mutations in mitochondria are neither advantageous nor disadvantageous (neutral mutations), but heteroplasmic cells incur a fitness cost proportional to the degree of heteroplasmy. Because no empirical data relate fitness to the degree of heteroplasmy, we consider three forms of fitness function to describe selection against heteroplasmy: concave, linear and convex (Fig 1A). For each fitness function, we vary the cost of heteroplasmy (c_h), given by c_h = 1 − h where h is the fitness of the most heteroplasmic cell in the population, to see how this affects the spread of U₁. We generate the concave fitness function by

the linear function by

and the convex function by

We also vary μ (mutation rate) and n (number of mitochondria) to ensure that our findings are robust. Second, we explore the effect of advantageous or deleterious mutations (non-neutral mutations) on the spread of U₁. Third, we relax the assumption of tight linkage between mating type and inheritance loci by exploring two cases: recombination between mating types and the absence of mating types altogether. Finally, we examine whether selection against heteroplasmy can explain the rare, but nevertheless important, exceptions to uniparental inheritance. To ensure that our results generalize to more than two mitochondrial types, we developed a second model that considers three mitochondrial types (S6 Model).

Results

When both mitochondrial haplotypes are neutral

We find that U₁ always replaces B₁, resulting in complete uniparental inheritance in the population (Fig 1B). These findings are independent of the number of mitochondria per cell (Fig 1C), mutation rate (Fig 1D), fitness function (Fig 1E), and cost of heteroplasmy (Fig 1F) (see S1–S10 Tables for more parameter combinations). We find the same results when we generalize the model to three mitochondrial haplotypes (S1 Fig).

General patterns

In our model, heteroplasmic cells are generated by mutation. During meiosis, heteroplasmic cells produce gametes with varying levels of heteroplasmy, including homoplasmic gametes. Uniparental inheritance maintains this variation created by meiosis, which leads to homoplasmic U₁B₂ cells (Fig 2A–2B and S2A–S2B Fig). Mutants that arise in U₁B₂ cells quickly segregate into U₁ gametes that carry mutant haplotypes only (Fig 3A–3B and S3A–S3B Fig), which leads to U₁B₂ cells that are homoplasmic for mutant mitochondria (Fig 2B and S2B Fig). Since we assume that mutations are neutral, cells homoplasmic for mutant mitochondria suffer no fitness costs.

U₁B₂ cells carrying mutant mitochondria produce B₂ gametes that also carry mutant mitochondria (Fig 3D and S3D Fig). When these B₂ gametes mate with B₁ gametes carrying wild type mitochondria, the resulting B₁B₂ cells are highly heteroplasmic (Fig 2C–2E and S2C Fig). As U₁ spreads, matings between U₁ and B₂ become more likely, increasing the level of heteroplasmy in both B₁B₂ cells and in B₁ and B₂ gametes (Figs. 2C–2E and 3C–3F and S2C and S3C–S3D Figs.). Increased levels of heteroplasmy reduce the fitness of both B₁ and B₂ gametes (w¯B1, w¯B2 in Fig 3A and S3A Fig) and B₁B₂ cells (w¯B1B2 in Fig 2A and S2A Fig). The difference in fitness between B₁ and B₂ becomes stronger (Fig 3A and S3A Fig) as more B₂ gametes that carry mutant mitochondria are produced (Fig 3D and S3D Fig). As a result U₁ spreads at the expense of B₁.

In the above description (Figs. 2 and 3), the mutation from B₁ to U₁ occurred in gametes homoplasmic for wild type mitochondria. When U₁ is introduced into heteroplasmic gametes, it takes fewer generations to reach equilibrium because B₂ gametes homoplasmic for mutant mitochondria are produced more quickly (S4 Fig). Our results are robust to changes in the frequency at which U₁ gametes are introduced (S5 Fig). For more detailed model dynamics, see S1 Text and S1–S2 Videos.

The effect of varying parameters

U₁ spreads more slowly when mutation rate (μ) is lower (Fig 1D) and number of mitochondria (n) is higher (Fig 1C). Reducing μ slows the spread of U₁ because mutant mitochondria are produced more slowly, slowing the generation of B₂ gametes that only carry the mutant haplotype. Increasing n has the same effect.

While varying the cost of heteroplasmy does not change the qualitative behavior of the model, it does affect the number of generations required for U₁ to replace B₁ (Fig 1F). In general, U₁ spreads more quickly when the cost of heteroplasmy is low for all three fitness functions (Fig 1F). Strong selection against heteroplasmy (e.g. c_h = 1) slows the production of B₂ gametes homoplasmic for the mutant haplotype because a transition via heteroplasmy is needed to lead to U₁B₂ cells homoplasmic for mutant mitochondria. Heteroplasmy levels thus remain low in B₁B₂ cells, and U₁ takes longer to replace B₁ (S6A and S6D Fig). At lower costs of heteroplasmy (e.g. c_h = 0.2), more B₂ gametes that are homoplasmic for the mutant haplotype are produced and levels of heteroplasmy in B₁B₂ cells increase, leading to a faster spread of U₁ (S6B and S6E Fig). Although levels of heteroplasmy in B₁B₂ cells increase even further as the cost of heteroplasmy approaches 0 (e.g. c_h = 0.01), selection against heteroplasmy is now very weak, which slows the spread of U₁ compared with c_h = 0.2 (S6C and S6F Fig). When the number of mitochondria is higher, U₁ spreads more quickly when the cost of heteroplasmy is low. This is because B₂ gametes homoplasmic for mutant mitochondria are produced more slowly at higher values of n and strong selection against heteroplasmy compounds this problem (S7 Fig). A similar logic can be applied to understand the differences between the three fitness functions. Since heteroplasmic cells are under weaker selection when fitness is concave (followed by linear and convex respectively) (Fig 1A), the level of heteroplasmy is highest using a concave function (S8 Fig). Thus, U₁ spreads more quickly using a concave function (followed by linear and convex respectively) when the cost of heteroplasmy is high because it is easier to generate heteroplasmic cells, and thus easier to generate B₂ gametes homoplasmic for mutant mitochondria, when selection against heteroplasmic cells is weaker (Fig 1F and S8 Fig). As the cost of heteroplasmy decreases, the number of generations for U₁ to spread under the three fitness functions converges because it becomes easier to generate B₂ gametes homoplasmic for mutant mitochondria (Fig 1F).

When mutations are deleterious

We next investigate how the U₁ allele spreads when mutations are non-neutral, as is the case for most mtDNA mutations [13]. We start by assuming that mutations are deleterious so that cells carrying mutant mitochondria are more strongly selected against than cells that carry wild type mitochondria. We assume that a mutation from wild type to mutant haplotype is more common than the reverse [14]. We let the probability of a mutation from mutant to wild type haplotype be μ_b = μ/100. We vary the selection coefficient of the mutant haplotype to see how this affects the spread of the U₁ allele (the fitness of a cell homoplasmic for the mutant haplotype is 1 − s_d, where s_d is the selection coefficient of the mutant haplotype). Essentially there are now two fitness functions: one governing the effect of mitochondria on cell fitness (where the selection coefficient determines the magnitude of the effect) and one governing the cost of heteroplasmy. For deleterious mutations, we assume that fitness decreases as a concave function of the number of mutants, as this relationship is experimentally established [15]. We examine both concave and convex fitness functions for selection against heteroplasmy (yielding two combinations).

Again, U₁ replaces B₁ unless the fitness of heteroplasmic cells and the fitness of deleterious mutants are governed by a concave function and the selection coefficient is sufficiently large (S9 Fig and S11–S12 Tables). U₁ generally spreads more slowly as s_d increases and it always spreads more slowly compared to when mutations are neutral (S11–S12 Tables). Stronger selection against mutant haplotypes leads to fewer B₂ gametes homoplasmic for mutant mitochondria, which slows the spread of U₁ (S10 Fig).

When mutations are advantageous

Next we explore the effect of advantageous mutations on the spread of U₁. In this case, cells that carry mutant haplotypes have an advantage over those carrying wild type haplotypes (the fitness of a cell homoplasmic for the wild type haplotype is 1 − s_a, where s_a is the selection coefficient of the mutant haplotype). We account for the rarity of advantageous mutations by setting μ_b = 100μ. Because it is unknown how fitness relates to the accumulation of advantageous mtDNA mutations, we model this relationship with both a concave and convex function. As in the deleterious case, we model selection against heteroplasmy by testing both concave and convex fitness functions (giving four combinations).

U₁ always replaces B₁ unless mutations are highly advantageous (s_a = 0.1) and both the fitness of heteroplasmic cells and the fitness of advantageous mutants are governed by a concave function (S9 Fig and S13–S14 Tables). U₁ spreads more quickly when s_a = 0.001 and s_a = 0.01 because B₂ gametes homoplasmic for mutant haplotypes now have a fitness advantage and are produced more quickly (S10 Fig). In contrast, U₁ spreads more slowly when s_a = 0.1 because the mutant haplotype quickly replaces the wild type as the dominant haplotype before U₁ has replaced B₁. Once B₁ gametes carry mostly mutant haplotypes, B₁ × B₂ matings are less costly because they predominantly involve mutant haplotypes. We find the same patterns for non-neutral mutations when we generalize our model to three mitochondrial types (S15 Table).

Recombination between mating type and inheritance loci

Previously, U × U matings were not possible because we assumed tight linkage between mating type and inheritance loci. But if we allow recombination to occur between these loci, U₁ × U₂ matings become possible. In this scenario, the number of gametes increases to four (B₁, B₂, U₁ and U₂), as does the number of genotypes (B₁B₂, U₁B₂, U₁U₂ and U₂B₁). There are three main ways in which mitochondrial inheritance could be regulated in U₁ × U₂ matings. (1) One U allele is dominant to the other, leading to uniparental inheritance; (2) each U allele ensures inheritance of its mitochondria, resulting in biparental inheritance; or (3) inheritance is more or less random so that some matings result in uniparental inheritance and some in biparental inheritance. We model all three cases.

When U₁ × U₂ matings lead to uniparental inheritance, the U₁U₂ genotype always spreads until it is fixed in the population, leading to complete uniparental inheritance (Fig 4A and S16–S18 Tables). When U₁ × U₂ matings lead to biparental inheritance, however, uniparental inheritance does not become fixed and the population reaches a polymorphic equilibrium (Fig 4B–4C). Under these conditions, the frequency of uniparental inheritance at equilibrium is ≤ 0.5 (S19–S21 Tables). Uniparental inheritance cannot exceed 0.5 because increasing the frequency of U₁ or U₂ simply increases the proportion of biparental U₁ × U₂ matings. The frequency of uniparental inheritance remains very low when we assume a concave fitness function (Fig 4B), but reaches its maximum (0.5) when we assume a linear or convex fitness function (Fig 4C) (see S12–S13 Figs. for an explanation).

When the probability of recombination (P_r) is sufficiently high (10⁻⁴ ≤ P_r ≤ 0.5 in S11 Fig), the U₁B₂ and U₂B₁ genotypes have the same frequency at equilibrium (S11B–S11D Fig). Now uniparental inheritance is no longer associated with a single mating type but is evenly split between the two mating types (S19–S21 Tables). When P_r is sufficiently small (P_r = 10⁻⁵ in S11 Fig), the recombination rate is so low that the mating type and inheritance loci are essentially linked and the U₁B₂ genotype becomes fixed (as in the general model) (S11A Fig).

When we assume a mixture of uniparental inheritance and biparental inheritance, we let U₁ × U₂ matings lead to biparental inheritance with probability P_b and to uniparental inheritance with probability 1 − P_b. Lowering P_b increases the frequency of uniparental inheritance, and uniparental inheritance becomes fixed when P_b = 0 (Fig 4A and 4E). Under linear and convex fitness functions, the equilibrium always maximizes the level of uniparental inheritance (Tables S22–S23). Under concave fitness, however, uniparental inheritance is only maximized for particular values of P_b (roughly P_b ≤ 0.2 for the parameter values we considered) (S22 Table; rows 2–3). (See S5 Model for how we determine when uniparental inheritance is maximized.)

We also find that uniparental inheritance can evolve in the complete absence of mating types. The no mating types scenario differs from the recombination case in that UB equals the sum of U₁B₂ and U₂B₁ at equilibrium (Fig 4A and 4F) (see S2 Text for more details).

Can selection against heteroplasmy explain the exceptions to uniparental inheritance?

In this section, we explore whether relaxing some of the assumptions in our general model can lead to mitochondrial inheritance patterns that resemble some of the known exceptions to uniparental inheritance. Exceptions to uniparental inheritance fall in three main categories: organisms that (1) regularly inherit mitochondria from both parents; (2) normally inherit mitochondria from one of the two parents but on occasion inherit mitochondria from both; and (3) inherit mitochondria from either or both parents.

Baker’s yeast, Saccharomyces cerevisiae, regularly inherits mitochondria from both parents (though uniparental inheritance also occurs), but heteroplasmy is transient because the diploid cell has only a few mitochondria [16] and divides repeatedly, which separates heteroplasmic cells into cells homoplasmic for either mitochondrial type (vegetative segregation) [9, 10]. Vegetative segregation is usually completed within twenty generations, but up to 50% of zygotes may be homoplasmic after the first division ([10] and references therein). Thus, Saccharomyces may restore homoplasmy as quickly as organisms that actively destroy one mitochondrial lineage [17]. Similarly, the geranium Pelargonium zonale often inherits cytoplasmic organelles from both parents (chloroplasts in this case). As with Saccharomyces, heteroplasmy is transient in Pelargonium because of rapid vegetative segregation of heteroplasmic cells shortly after syngamy [9]. We added mitotic divisions to our model to test whether vegetative segregation could maintain biparental inheritance under selection against heteroplasmy. When we include mitosis before selection (which assumes that vegetative segregation occurs swiftly, before selection has time to act), uniparental inheritance does not spread, provided that the number of mitochondria is low (n = 4) and the number of divisions is high (S24 Table; rows 7–8). Under these conditions, biparental inheritance is stable because heteroplasmic cells resulting from biparental inheritance segregate into homoplasmic cells before selection acts. If there are insufficient mitotic divisions, or if selection acts before vegetative segregation is complete, then uniparental inheritance replaces biparental inheritance, although it spreads much more slowly than when there are no mitotic divisions (S24 (rows 3–6) and S25 Tables). When there are more mitochondria per cell (e.g. n = 8), biparental inheritance is only stable if the number of cell divisions increases to compensate (S24 Tables; rows 9–10). Thus, biparental inheritance can be stable under selection against heteroplasmy but only under a narrow set of conditions, explaining why this form of inheritance is rare.

In other isogamous organisms, including the acellular slime molds Physarum polycephalum and Didymium iridis and the algae Chlamydomonas reinhardtii, mitochondria from both gametes mix before one mitochondrial lineage is destroyed post-fertilization, often by nucleases [18–20]. This mechanism is not perfect and these organisms sometimes deviate from strict uniparental inheritance [9, 18–20]. While uniparental inheritance is the norm in the slime mold P. polycephalum, sometimes both mitochondrial lineages survive, leading to varying degrees of biparental inheritance [18]. Could uniparental inheritance still spread under such conditions? Since mating types and inheritance loci are tightly linked in Physarum [18], we explore this question using our general model that assumes linkage. Now, U₁ × B₂ matings lead to biparental inheritance with probability P_b and to uniparental inheritance with probability 1 − P_b. For the parameter values that we examined, the U₁B₂ genotype always goes to fixation when P_b < 1 and the fitness function is linear or convex (S26 Table). (When fitness is concave, P_b must be roughly <0.05 for the U₁B₂ genotype to become fixed.) Under these conditions, the frequency of biparental inheritance at equilibrium is equal to P_b (S26 Table). In this scenario, the level of biparental inheritance in the population simply reflects the likelihood that an individual mating results in biparental inheritance.

Chlamydomonas reinhardtii and Didymium iridis can inherit mitochondria from either or both parents [19, 20]. Chlamydomonas normally inherits mitochondria from the mt − parent and chloroplasts from the mt + parent, but under some circumstances it can inherit mitochondria from mt + and chloroplasts from mt − or mitochondria and chloroplasts from both [20]. Didymium iridis has random, biased, or dominant patterns of uniparental inheritance. Under random uniparental inheritance, either parental strain is equally likely to be the mitochondrial donor while, under biased inheritance, one strain is more likely to be the mitochondrial donor [19]. Under dominant inheritance, one strain is always the donor. Didymium also has low levels of biparental inheritance [19]. In this scenario, we test whether selection against heteroplasmy could lead to the evolution of a system with a mixture of uniparental inheritance (from either parent) and biparental inheritance. We assume that mating types can recombine and that U₁ × U₂ matings can lead to mitochondria being inherited from U₁, U₂ or both. Mitochondria are inherited from U₁ with probability PU1, from U₂ with probability PU2 and from both parents with probability P_b (where PU1+PU2+Pb=1). Now, uniparental inheritance comes from U₁ × B₂ matings, U₂ × B₁ matings and those U₁×U₂ matings with uniparental inheritance. Irrespective of the values of PU1 and PU2, we find the same results as with our earlier model in which U₁×U₂ matings led to a mixture of uniparental and biparental inheritance (S22–S23 Tables). This is because equilibrium depends only on the value of P_b. (Since uniparental inheritance quickly eliminates most heteroplasmic cells, U₁U₂ cells are almost entirely homoplasmic regardless of which gamete donates mitochondria.) Consequently, different probabilities of inheriting mitochondria biparentally (P_b), from mating type 1 (PU1) or from mating type 2 (PU2) lead to a range of inheritance patterns that include uniparental inheritance (from both parents) and biparental inheritance (see S27 Table for some examples).

Lastly, selection against heteroplasmy provides an explanation for the cases in which mitochondria are inherited from one parent while chloroplasts are inherited from the other (e.g. in Chlamydomonas and pines [20, 21]). If uniparental inheritance simply evolved to maintain homoplasmy in cells, it should not matter which parent donates mitochondria or chloroplasts.

Discussion

Our model shows that selection against heteroplasmy can lead to the fixation of uniparental inheritance in an ancestrally biparental population. We find that uniparental inheritance replaces biparental inheritance under almost all tested scenarios and parameter values. Our model also explains many of the known exceptions to strict uniparental inheritance. We show that uniparental inheritance can replace biparental inheritance whether mutations lead to neutral or non-neutral haplotypes. Relaxing our initial assumptions of pre-existing mating types and lack of recombination does not prevent uniparental inheritance from evolving. As we make no attempt to resolve the evolution of mating types within the context of mitochondrial inheritance, as others have previously attempted [1, 22], our findings thus leave open the possibility that mating types preceded uniparental inheritance, evolved as a consequence of uniparental inheritance, or evolved after uniparental inheritance.

In contrast to previous models, we show that uniparental inheritance can spread under realistic mutation rates and number of mitochondria per cell. The lowest value of μ that we tested (10⁻¹⁰) is eight orders of magnitude lower than required by the genomic conflict theory [1] and compares favorably with empirical mutation rates (10⁻⁷ to 10⁻⁸ per site per generation [23–25]). Both the genomic conflict and mutation clearance hypotheses require unrealistic mutation rates and number of mitochondria per cell for uniparental inheritance to replace biparental inheritance, while uniparental inheritance cannot replace biparental inheritance under any parameter values in the mitochondrial-nuclear coadaptation model [1]. The genomic conflict model requires a mutation rate of 1% per generation before uniparental inheritance can replace biparental inheritance [1]. The only known example that satisfies this assumption is the petite mutant in Saccharomyces cerevisiae, which is a hyper-mutable selfish mitochondrion that can spontaneously arise at a rate of 1% per generation [26]. Under this mutation rate, however, the genomic conflict model requires that cells contain at least 50 mitochondria [1], whereas most extant isogamous species, including Saccharomyces, contain fewer than 20 mitochondria at syngamy [16, 18]. As mutant mitochondria lack a transmission advantage over wild type mitochondria in the mutation clearance hypothesis, the mutation clearance model requires even higher mutation rates [1]. To the best of our knowledge, no extant organism satisfies the assumptions of the genomic conflict or mutation clearance hypotheses.

Why do our results differ from the findings of previous models? In the genomic conflict and mutation clearance models, wild type mitochondria mutate to selfish or deleterious mitochondria. Biparental inheritance results in cells that are heteroplasmic for wild type and mutant mtDNA, while U₁ gametes mostly contain wild type mitochondria [1]. Because U₁ purges B₂ gametes of mutant mitochondria, B₁×B₂ matings involve increasingly fewer mutant mitochondria as the frequency of U₁ increases [1, 5]. U₁ is thus subject to negative frequency-dependent selection, and the population reaches equilibrium well before uniparental inheritance replaces biparental inheritance at realistic mutation rates [1]. The mitochondrial-nuclear coadaptation model assumes that mitochondria are well matched or poorly matched to nuclear alleles [1, 8]. Because mutation can lead to matched nuclear-mitochondrial states becoming unmatched, the effective mitochondrial mutation rate is lower in the mitochondrial-nuclear coadaptation model, which prevents uniparental inheritance from displacing biparental inheritance under any parameter values [1].

Evidence for a cost of heteroplasmy comes from a recent study that compared the effect of two mtDNA haplotypes (NZB and 129S6) in a cogenic nuclear background on the functioning of mice [12]. Mice homoplasmic for NZB or 129S6 were phenotypically normal, but NZB-129S6 heteroplasmic mice suffered from reduced activity, lowered food intake, compromised respiration, heightened stress response, and impaired cognition [12]. While the mechanism(s) behind the cost of heteroplasmy is unknown, there are a few possibilities. Heteroplasmy may disrupt cell signaling by altering production of reactive oxygen species (ROS) [27] and there are indications that heteroplasmy can increase mitochondrial ROS levels [28, 29], leading to phenotypes that differ from cells that are homoplasmic for either haplotype [29, 30]. Alternatively heteroplasmy may lead to deleterious interactions between polypeptides from different mitochondria within the same electron transport chain [12, 31]. Because chloroplasts also contain independent genomes, are involved in cellular bioenergetics, and generally show uniparental inheritance [9], our findings likely apply to both mitochondria and chloroplasts.

Although the evidence in mice is compelling [12], it is unknown whether selection against heteroplasmy is a general phenomenon in eukaryotes. While Sharpley and colleagues [12] used different mitochondrial lineages to construct heteroplasmic individuals, our model assumes that mutations accumulated within a single generation can cause mitochondrial types to become sufficiently distinct to lead to negative effects for the cell. At this stage we do not know how different mitochondrial genomes have to be for selection against heteroplasmy to apply. It could also be that there are regions of the genome in which heteroplasmic mutations have a stronger effect on fitness than others. To support or refute our model, we now need solid empirical data on a range of organisms showing the cost, if any, of heteroplasmy on organism fitness.

While we have referred to n as the number of mitochondria in the cell, n actually refers to the number of segregating units of mtDNA at syngamy. Mitochondria pack DNA into DNA-protein complexes called nucleoids, which themselves may contain multiple copies of mtDNA [32, 33]. It is currently unknown whether the segregating unit is the mtDNA molecule itself, the nucleoid, the mitochondrion or another level of mtDNA organization [33]. But as nucleoids are predominantly homoplasmic, even in heteroplasmic tissues [32, 33], the number of mitochondria may be a reasonable approximation of the number of segregating units in the cell. If the segregating unit is at a lower level of organization (e.g. the mtDNA molecule), then n, as used in our model, will apply to the number of segregating units not the number of mitochondria per cell (e.g. n = 200 would then apply to a cell with 200 segregating units, which may be a cell with far fewer than 200 mitochondria).

By assuming an infinite population size, a common assumption in studies of this kind [1, 5, 6, 8] we have ignored genetic drift, which can be a powerful force in population genetics. While it is beyond the scope of this study to formally model the effects of genetic drift on the evolution of uniparental inheritance, we can anticipate some of its effects. As the mutation leading to uniparental inheritance has a small advantage when its frequency is low, genetic drift will lead to the frequent loss of those mutations. Thus, the initial invasion of a mutation for uniparental inheritance may be largely determined by genetic drift rather than by positive selection. As the frequency of uniparental inheritance increases, however, so too does its advantage, reducing the probability that the mutation is lost to drift. The potential for rare mutations to be lost to drift is not unique to our model. The genomic conflict hypothesis requires stringent conditions for uniparental inheritance mutations to invade [6, 34]. Under this hypothesis, a mutation for uniparental inheritance must arise within a population that contains selfish mutants but in which the selfish mutant is not fixed. Otherwise, uniparental inheritance cannot become associated with non-selfish mitochondria. Any mutations leading to uniparental inheritance that arise outside of this window will have no selective advantage and will be more likely to be lost by genetic drift [6, 34].

Conclusion

Selection against heteroplasmy has implications for the evolution of the mitochondrial genome. Because of a smaller effective population size, which is more strongly affected by genetic drift, and higher mutation rates, mtDNA should be less conserved than the nuclear genome [35, 36]. Indeed, mitochondrial transfer RNAs and synonymous sites mutate 5–50 times more frequently than comparable elements in the nuclear genome [35, 37]. Because the mitochondrial genome is effectively asexual, any deleterious mutations in the fittest haplotype cannot be rescued (except by unlikely back mutations). This effect, known as Muller's Ratchet, should eventually lead to irreparable genome meltdown [38, 39]. In stark contrast to theoretical predictions, however, mitochondrial coding genes are more conserved than analogous nuclear oxidative phosphorylation genes [36]. When mtDNA mutates, only one of the many mtDNA molecules in the cell is affected, leading to a heteroplasmic cell. Selection against heteroplasmy should reduce the probability that mtDNA mutations spread throughout the cell, which, in turn, should oppose changes to mtDNA. Thus, selection against heteroplasmy may not only explain the evolution of uniparental inheritance but also why mitochondrial coding genes have thus far managed to resist the effects of Muller's Ratchet.

Model

Our model tracks the distribution of cell types through each stage of the life cycle across multiple generations. The redistribution of cell types is based on probability theory, but the model itself is deterministic. We assume that the population is effectively infinite and unaffected by genetic drift, as is regularly assumed in models such as ours [1, 5, 6, 8]. Consequently, the probability that a cell takes a particular state equates to the proportion of that cell type in the population. We take a similar approach to previous models [1, 5], but our model differs slightly in our treatment of mutation. Hastings does not include mutation [5], while Hadjivasiliou and colleagues treat mutation as a one-way process from wild-type to mutant mitochondria in the conflict and mutation clearance models [1]. When examining the mitochondrial-nuclear coadaptation model, however, Hadjivasiliou and colleagues allow mutation to proceed both ways as we have done here [1]. In our model, mutation is designed to capture the ability of a mitochondrial type to mutate from its current state to other haplotypes (one type in our main model and two types in our supplementary model, but an extremely large number of haplotypes in reality).

Diploid cell types are described by the vector Mt,τα=(i,G), where i corresponds to the number of mutant mitochondria and takes values in {0,1…n}, t indicates the generation, and τ_α indicates the stage of the life cycle. If we know the number of mutant mitochondria (i), the number of wild type mitochondria (which we denote j) is fixed as j = n – i. G indicates the nuclear genotype and takes values in {U₁B₂,B₁B₂}. Gametes are described by the vector Mt,τα=(p,g), where p is the number of mutant mitochondria and takes values in {0,1…n / 2} and g represents the nuclear allele and takes values in {U₁,B₁,B₂}. The probability of obtaining a particular diploid cell type is written as P(Mt,τα=(i,G)) and the probability of obtaining a particular gamete is written as P(Mt,τα=(p,g)).

These probabilities can also be thought of as the proportion of the population with that particular cell or gamete type.

There are n+1 total mitochondrial states for diploid cells and n / 2+1 possible mitochondrial states for haploid cells. For the case in which mating type and inheritance loci are linked, the total number of diploid cell types is 2(n+1) while the total number of haploid cell types is 3(n / 2+1). We obtained numerical solutions to our model via scripts that we developed in MATLAB (version 2013b).

Initialization

The starting population is evenly split between B₁ and B₂ gametes, and all gametes contain type wild type mitochondria (i.e. P(M0,τ1=(0,B1))=0.5,P(M0,τ1=(0,B2))=0.5 and P(M0,τ1=(p,g))=0,∀p>0andg=U1). We first allow this population to reach equilibrium, which we define as the point at which the proportion of cell types change by less than 10⁻¹² (except when the probability that a mitochondrion mutates into another mitochondrion is 10⁻¹⁰ (μ = 10⁻¹⁰), in which case we define equilibrium to be a change of less than 10⁻¹³). We then introduce U₁ gametes that are homoplasmic for wild type mitochondria by setting P(Mge1,τ1=(0,U1))=0.01, where ge1 is the number of generations taken to reach the first equilibrium. To maintain the total proportion of the population at 1, we remove the corresponding proportion of cells from the B₁ population (i.e.P(Mge1,τ1=(0,U1))=P(Mge1,τ1=(0,U1))−0.01). In two instances, we alter the way in which U₁ is introduced. In S4 Fig, we introduce U₁ into the most heteroplasmic gamete with a frequency greater than 0.01, and in S5 Fig we vary the introductory frequency of U₁. Our life cycle is very similar to the life cycle used by Hadjivasiliou and colleagues [1], which examined the genomic conflict, mutational clearance, and mitochondrial-nuclear coadaptation hypotheses.

Random mating

Gametes with n / 2 mitochondria randomly mate with the opposite mating type to produce diploid cells containing n mitochondria. In effect, this is random mating in which all matings between the same mating type (i.e. U₁U₁, B₁B₁, B₂B₂ and U₁B₁) are lethal, and the only viable genotypes are U₁B₂ and B₁B₂.

Biparental mating

Consider a biparental mating involving a gamete in state Mt,τ1=(p,B1), where τ₁ is the gamete stage of the life cycle. For this gamete to produce a diploid cell with type Mt,τ2=(i,B1B2), where τ₂ is the diploid cell stage of the life cycle that precedes mutation, it must mate with a gamete of type Mt,τ1=(i−p,B2). The probability of this mating is 2P(Mt,τ1=(p,B1))P(Mt,τ1=(i−p,B2)), where the factor of 2 accounts for the two ways in which we can choose B₁ and B₂ (B₁ then B₂ or B₂ then B₁). We restrict the values of p and i – p to biologically valid combinations. First, 0 ≤ p ≤ n / 2, as the B₁ gamete cannot carry negative numbers of mutant mitochondria nor can it contain more mutant mitochondria than the total number of mitochondria in the gamete. Likewise, 0 ≤ i – p ≤ n / 2 for the B₂ gamete, which, when rearranged, gives i – (n / 2) ≤ p ≤ i. Valid values for p lie in the range of intersection of these two inequalities, giving max(0,i – (n / 2)) ≤ p ≤ min(n / 2,i).

We can thus obtain the probability of forming any given diploid cell type after random mating with the sum,

Uniparental mating

Because uniparental matings between U₁ and B₂ gametes contain mitochondria from U₁ alone, U₁B₂ cells initially have n / 2 mitochondria. To restore the total complement of n mitochondria, we sample n / 2 mitochondria with replacement from the n / 2 mitochondria in the U₁B₂ cell and add the n / 2 sampled mitochondria to the original set of mitochondria to form a cell with n mitochondria.

For a gamete with identity Mt,τ1=(p,U1) to produce a diploid cell with identity Mt,τ2=(i,U1B2), it must sample n / 2 mitochondria containing i – p mutant mitochondria and n / 2 – (i – p) wild type mitochondria. The mitochondrial state of the B₂ gamete is irrelevant because its mitochondria are discarded and we will refer to this cell as Mt,τ1=(r,B2).

Sampling of the n / 2 mitochondria follows a binomial distribution, which we denote T(i – p;n / 2,(2p) / n), where i – p refers to the number of mutant mitochondria that need to be sampled, n / 2 refers to the number of mitochondria being sampled, and (2p) / n is the probability of drawing a single mutant mitochondrion from a U₁B₂ cell with p (out of n / 2) mutant mitochondria (where (2p) / n is obtained by rearranging p / (n / 2)).

The probability of sampling i – p mutant mitochondria (and (n / 2) – (i – p) wild type mitochondria) is given by

The restrictions on p and i – p are the same as those in biparental mating. Since U₁ will form the same initial U₁B₂ cell regardless of the B₂ gamete with which it mates, the probability of producing each type of U₁ gamete is multiplied by the probability of selecting any B₂ gamete. The probability of forming a given U₁B₂ cell after random mating is determined by

Mutation

We denote the post-mutation states of cells as Mt,τ3=(i,G), (where τ₃ indicates the post-mutation life cycle stage). If we define the number of wild type mitochondria that mutate to mutant mitochondria to be a and the number of mutant mitochondria that mutate to wild type mitochondria as b, a post-mutation cell in state Mt,τ3=(i,G) must be derived from a pre-mutation cell in state Mt,τ2=(i−a+b,G) (because the pre-mutation cell gains a mutant mitochondria and loses b mutant mitochondria to form the post-mutation cell). Similarly, if the post-mutation cell has j wild type mitochondria, then the pre-mutation cell must have j + a – b wild type mitochondria, where j = n – i.

First, we must work out the probability that a cell mutates a of its wild type mitochondria to mutant mitochondria. We define Y(a;n – i + a – b,μ) as the probability that a pre-mutation cell has a mutations in its n – i + a – b wild type mitochondria given that each mitochondrion mutates with probability μ. The accumulation of mutations is binomially distributed such that

Similarly, we define Y(b;i – a + b,μ_b) to be the probability that a pre-mutation cell acquires b mutations in its i – a + b mutant mitochondria given that each mitochondrion mutates with probability μ_b. This probability is given by

For any combination of values for a, b and i, multiplying Y(a;n – i + a – b,μ) by Y(b;i – a + b,μ_b) gives the probability of a particular transition from a pre-mutation cell with identity Mt,τ3=(i−a+b,G) to a post-mutation cell with identity Mt,τ3=(i,G). To get the overall probability that such a transition occurs, we multiply the probability of the transition by the proportion of pre-mutation cells in the population. To produce the post-mutation population, we sum all possible transitions between pre-mutation and post-mutation cells. All valid transitions must satisfy 0 ≤ a ≤ i (because the post-mutation cell cannot receive more than i mutant mitochondria) and 0 ≤ b ≤ n – i (because the post-mutation cell cannot receive more than n – i wild type mitochondria). Thus, we can determine the post-mutation population by

In the neutral scenario, μ = μ_b (i.e. the rate of mutation from wild type to mutant is equal to the rate of mutation from mutant to wild type).

Selection

The relative fitness of a cell, w(i), is a measure of how likely a cell type is to survive and reproduce, and we assume that cells carrying multiple mitochondrial types have lower fitness. For the first fitness function, the relative fitness of a cell with i mutant mitochondria is determined according to the following piecewise concave function:

for even values of n and 0 ≤ c_h ≤ 1, where c_h is the cost of heteroplasmy. In this function, a cell containing n / 2 mutant and n / 2 wild type mitochondria has minimum relative fitness.

The post-selection population of each cell type is then given by:

We also make use of two alternative fitness functions. The first of these is the piecewise linear function:

The third fitness function is the piecewise convex function:

We normalize the post-selection population by

where

so that the sum of the proportions of the population equals 1.

Meiosis

The cell must first duplicate its chromosomes and double its mitochondrial complement (from n to 2n). This cell with 2n mitochondria then produces gametes with n / 2 mitochondria. Meiosis occurs in two steps. First, we sample n mitochondria with replacement from a cell containing n mitochondria and add the set of sampled mitochondria to the original set of mitochondria to form a cell containing 2n mitochondria (this is the same process that occurs in uniparental mating only with n mitochondria rather than n / 2 mitochondria). We let Mt,τ6=(l,2G) represent the cell with doubled mitochondria and nuclear genotype, where l takes values in {0,1…2n} and 2G takes values in {U₁U₁B₂B₂,B₁B₁B₂B₂}.

For a cell to contain l mutant mitochondria after duplication of its mitochondria, it must sample l – i mutant mitochondria. We denote the probability of sampling l – i mutant mitochondria from Mt,τ5=(i,G) as F(l – i;n,i / n). Sampling follows a binomial distribution such that

We obtain Mt,τ6=(l,2G) by

During the second step of meiosis, the cells with 2n mitochondria produce gametes with n / 2 mitochondria. Biologically, this occurs in two steps. In meiosis 1, the homologous chromosomes are pulled apart to produce two haploid cells that contain two identical nuclear alleles (sister chromatids) and n mitochondria. In meiosis 2, the two cells divide to produce four gametes, each with a single nuclear allele and n / 2 mitochondria. Since mitochondria segregate independently of nuclear alleles during cell partitioning, we model this as a single step.

We define S(p;2n,l,n / 2) to be the probability of obtaining p mutant mitochondria in n / 2 draws from a cell in state Mt,τ6=(l,m,2G). Here, sampling is without replacement and follows a hypergeometric distribution, giving

Gametes produced by meiosis are represented by Mt+1,τ1=(p,g). We determine the probability of obtaining a particular gamete using

and

Factors of 1 / 2 in the above three equations take into account that half of the gametes produced from parent cells with nuclear genotype U₁B₂ will carry the U₁ allele and the other half will carry the B₂ allele (with similar logic applied for gametes produced from parent cells with nuclear genotype B₁B₂). Meiosis completes a single generation of the life cycle.

Relative fitness of cells

The relative fitness of U₁B₂ cells is given by

while the relative fitness of B₁B₂ cells is

Relative fitness of gametes

Although gametes are not subject to selection in our model, and thus do not technically have fitness values, it is informative to track gamete relative fitness throughout the simulation. We define a gamete’s relative fitness as the fitness that a diploid cell would have if it had the same mitochondrial composition as the gamete. Since gametes contain n / 2 mitochondria, they will have minimum fitness when they carry n / 4 wild type and n / 4 mutant mitochondria. To rescale the fitness function, we substitute n / 2 for n in the diploid cell fitness functions. For example, Equation (1) becomes

Once the fitness function is scaled to gametes, we can determine the relative fitness of the three gametes by

and

See S1–S6 Model for details of the other models.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14

Attachment 15

Attachment 16

Attachment 17

Attachment 18

Attachment 19

Attachment 20

Attachment 21

Attachment 22

Attachment 23

Attachment 24

Attachment 25

Attachment 26

Attachment 27

Attachment 28

Attachment 29

Attachment 30

Attachment 31

Attachment 32

Attachment 33

Attachment 34

Attachment 35

Attachment 36

Attachment 37

Attachment 38

Attachment 39

Attachment 40

Attachment 41

Attachment 42

Attachment 43

Attachment 44

Attachment 45

Attachment 46

Attachment 47

Attachment 48

Attachment 49

Attachment 50

Attachment 51

Attachment 52

Attachment 53

Attachment 54

Attachment 55

Attachment 56