We will consider three different classes of proposed generating mechanisms for the mtDNA bottleneck: those proposed in Cao et al. (2007); Cree et al. (2008) and Wai et al. (2008). We will refer to these mechanisms by their leading author name. The Cree mechanism involves random replication and degradation of mtDNAs throughout development, and binomial partitioning of mtDNAs at cell divisions. The Cao mechanism involves partitioning of clusters of mtDNA at each cell division, thus providing strong stochastic effects associated with each division. We consider a general set of dynamics through which this cluster inheritance may be manifest, including the possibility of heteroplasmic ‘nucleoids’ of constant internal structure (Jacobs et al., 2000), sets of molecules or nucleoids within an organelle, homoplasmic clusters, and different possible cluster sizes (see Appendix 1). The Wai mechanism involves the replication of a subset of mtDNAs during folliculogenesis. We note that this latter mechanism can be manifest in several ways: (a) through slow random replication of mtDNAs (so that, in any given time window, only a subset of mtDNAs will be actively replicating) or (b) through the restriction of replication to a specific subset of mtDNAs at some point during development. We will refer to these different manifestations as Wai (a) and Wai (b) respectively. The Wai (a) mechanism and the Cree model can both be addressed in the same mechanistic framework (with potentially different parameterisations): if the rate of random replication in the Cree model is sufficiently low during folliculogenesis, only a subset of mtDNAs will be actively replicating at any given time during this period, thus recapitulating the Wai (a) mechanism (see Appendix 1). We will henceforth combine discussion of the Wai (a) and Cree mechanisms into what we term the birth-death-partition (BDP) mechanism.

We seek a physically motivated mathematical model for the bottleneck that is capable of reproducing each of these mechanisms. Our general model for the bottleneck (detailed description in ‘Materials and methods’) involves a ‘bottom-up’ representation of mtDNAs as individual intracellular elements capable of replication and degradation (Figure 1B) with rates λ and ν respectively. A parameter S determines whether these processes are deterministic (specific rates of proliferation) or stochastic (replication and degradation of each mtDNA is a random event). These rates of replication and degradation of mtDNA are likely strongly linked to mitochondrial dynamics within cells, through the action of mitochondrial quality control (Twig et al., 2008; Hill et al., 2012) modulated by mitochondrial fission and fusion (Detmer and Chan, 2007; Youle and van der Bliek, 2012; Hoitzing et al., 2015), which can act to recycle weakly-perfoming mitochondria (Mouli et al., 2009; Twig and Shirihai, 2011). This quality control can be represented through the degradation rates assigned to each mtDNA species, which may differ (for selective quality control) or be identical (for non-selective turnover).

The proportion of mtDNAs capable of replication is controlled by a parameter α in our model, dictating the proportion of mtDNAs that may replicate after a cutoff time T. Thus, if α = 1, all mtDNAs may replicate; if α < 1, replication of a subset proportion α of mtDNAs is enforced at this cutoff time. At cell divisions, mtDNAs may be partitioned either deterministically, binomially, or in clusters according to a parameter c (Figure 1C).

The copy number of mtDNA per cell is observed to vary dramatically during development, with dynamic phases of copy number depletion and different rates of subsequent recovery observed. Additionally, cell divisions occur in the germline at different rates during development, with cells becoming largely quiescent after primary oocytes develop. To explicitly model these different dynamic regimes, and the behaviour of mtDNA copy number during each, we include six different dynamic phases throughout development, each with different rates of replication and degradation (labelled with subscript i labelling the dynamic phase: hence λ₁, ν₁,…,λ₆, ν₆), and allowing for different rates of cell division or quiescence. This protocol enables us to explicitly model effects of changing population size throughout development rather than assuming dependence on a single, coarse-grained effective population size; and to include the effects of specific and varying cell doubling times. A summary of symbols used in our model and throughout this article is presented in Figure 1D.

Our model, with suitable parameterisation, can thus mirror the dynamics of the Cree and Wai (a) mechanisms (stochastic dynamics and binomial partitioning, which we refer to as the BDP mechanism); the Cao mechanism (clustered partitioning); and Wai (b) mechanism (deterministic dynamics, restricted subset of replicating mtDNAs). The Cao mechanism, partitioning of clusters of mtDNA molecules, represents the expected case if mtDNA is partitioned in colocalised ‘nucleoids’ within each organelle (or in other sub-organellar groupings). The size of mtDNA nucleoids is debated in the literature (Bogenhagen, 2012; Kukat and Larsson, 2013; Wallace and Chalkia, 2013) (although recent evidence from high-resolution microscopy suggests that nucleoid size is generally <2 (Jakobs and Wurm, 2014), consonant with recent evidence that individual nucleoids may be homoplasmic [Poe et al., 2010]); our model allows for inheritance of homoplasmic or heteroplasmic nucleoids of arbitrary characteristic size c, thus allowing for a range of sub-organellar mtDNA structure. We discuss the impact of mixed or fixed nucleoid content in Appendix 1.

We take data on mtDNA copy number in germ line cells in mice from three experimental studies (Cao et al., 2007; Cree et al., 2008; Wai et al., 2008). We also use data from two experimental studies on heteroplasmy variance in the mouse germ line during development (Jenuth et al., 1996; Wai et al., 2008). These heteroplasmy variance studies employ intracellular combinations of the same pairing of mtDNA haplotypes (NZB and BALB/c), modelling two different mtDNA types within a cellular population. These data, by convention (Samuels et al., 2010), are normalised by heteroplasmy level h, giving

where normalised variance $\mathbb{V}\prime \left(h\right)$ is a quantity that will be often used subsequently. This normalised variance controls for the effect of different or changing mean heteroplasmy, and thus allows a comparison of heteroplasmy variance among samples with different mean heteroplasmies and subject to heteroplasmy change with time. We use a time of 100 dpc to correspond to mature oocytes (see ‘Materials and methods’). We take data on cell doubling times from a classical study (Lawson and Hage, 1994) (see ‘Materials and methods’). A possible summary of these data (although they provoke ongoing debate; see ‘Discussion’) is that, as shown in Figure 2A, the existing data on normalised heteroplasmy variance shows initially low variance until ∼7.5 dpc (days post conception, which we use as a unit of time throughout), rising to intermediate values between 7.5 and 21 dpc, gradually rising further subsequently to become large in the mature oocytes of the next generation. In Figure 2A, and throughout this article, experimentally measured data will be depicted as circular or polygonal points, and inferred theoretical behaviour will be depicted as lines or shaded regions.

Figure 2

Download asset Open asset

Figure 2A shows mtDNA population dynamic trajectories resulting from optimised parameterisations of each of the mechanisms we consider (see ‘Materials and methods’). In Figure 2B we show posterior probabilities on each of these mechanisms. These posterior probabilities give the inferred statistical support for each mechanism, derived from model selection performed with approximate Bayesian computation (ABC) (Beaumont et al., 2002; Toni et al., 2009; Sunnåker et al., 2013; Johnston, 2014) using uniform priors. ABC involves choosing a threshold value dictating how close a fit to experimental data is required to accept a particular model parameterisation as reasonable. In our case, this goodness-of-fit is computed using a comparison of squared residuals associated with the trajectories of mean mtDNA copy number and normalised heteroplasmy variance (see ‘Materials and methods’ and Appendix 1). Each of the experimental measurements corresponds to a sample variance, derived from a finite number of samples of an underlying distribution of heteroplasmies, and therefore has an associated uncertainty and sampling error (Wonnapinij et al., 2010). The reasonably small sample sizes used in these sample variance measurements are likely to underestimate the underlying heteroplasmy variance (the target of our inference). Our ABC approach naturally addresses these uncertainties by using summary statistics derived from sampling a set of stochastic incarnations of a given model, where the size of this set is equal to the number of measurements contributing to the experimentally-determined statistic (see ‘Materials and methods’). Figure 2B clearly shows that as the ABC threshold is decreased, requiring closer agreement between the distributions of simulated and experimental data, the posterior probability of the BDP model increases, to dramatically exceed those of the other models. This increase indicates that the BDP model is the most statistically supported, and capable of providing the best explanation of experimental data (which can be inutitively seen from the trajectories in Figure 2A). We note that ABC model selection automatically takes model complexity into account, and conclude that the BDP mechanism is the best supported proposed mechanism for the bottleneck. Briefly, this result arises because the BDP model produces increasing variance both due to early cell division stochasticity and later random turnover. By contrast, the Cao model alone only increases variance in early development when cell divisions are occurring. Qualitatively, this behaviour through time holds regardless of cluster (nucleoid) size and regardless of whether clusters are heteroplasmic or homoplasmic (allowing heteroplasmic clusters decreases the magnitude of heteroplasmy variance but not its behaviour through time, see Appendix 1). The Wai (b) model alone similarly only increases variance at a single time point (later, during folliculogenesis).

In Wai et al. (2008), visualisations of cells after BrU incorporation show that a subset of mitochondria retain BrU labelling, which the authors suggest indicates that a subset of mtDNAs are replicating. In Appendix 1, we show that the BDP model also results in the observation of only a subset of replicating mtDNAs over the timeframe corresponding to these experimental results. These observations thus correspond to results expected from the random turnover from the BDP model. We also note the mathematical observation that the Wai (b) mechanism requires the replication of <1% of mtDNAs during folliculogenesis to yield reasonable heteroplasmy variance increases (Figure 2A shows the optimal case with α = 0.006; optimal fits to data generally show 0.005 < α < 0.01), and the proportions of loci visible in Wai et al. (2008) are substantially higher than this required 1% value.

We show in Appendix 1 that the heteroplasmy statistics corresponding to binomial partitioning also describe the case where the elements of inheritance are heteroplasmic clusters, where the mtDNA content of each cluster is randomly sampled from the population of the cell (either once, as an initial step, or repeatedly at each division). This similarity holds broadly, regardless of whether the internal structure of clusters is constant across cell divisions or allowed to mix between divisions. The BDP model, in addition to describing the partitioning of individual mtDNAs, also thus represents the statistics of mtDNA populations in which heteroplasmic nucleoids are inherited (Jacobs et al., 2000), or individual organelles containing a mixed set of mtDNAs or nucleoids are inherited, regardless of the size of these nucleoids (see ‘Discussion’).

Having used ABC model selection to identify the BDP model as the most statistically supported, we can also use ABC to infer the values of the governing parameters of this model given experimental data. Figure 3A,B shows the trajectories of mean copy number and mean heteroplasmy variance resulting from model parameterisations identified through this process. Figure 3C shows the inferred behaviour of mtDNA degradation rate ν in the model, a proxy for mtDNA turnover (as the copy number is constrained). Turnover is generally low during cell divisions, allowing heteroplasmy variance to increase due to stochastic partitioning. Turnover then increases later in germ line development, resulting in a gradual increase of heteroplasmy variance after birth until the mature oocytes form in the next generation.

Figure 3

Download asset Open asset

Figure 3D shows posterior distributions on the copy number minimum and total turnover (see ‘Materials and methods’) resulting from this process; posteriors on all other parameters are shown in Appendix 1. Substantial flexibility exists in the magnitude of the copy number minimum, illustrating that observed heteroplasmy variance can result from a range of bottleneck sizes from ∼200 to >10³; going some way towards reconciling the conflict between Cao et al. (2007) and Cao et al. (2009) and Cree et al. (2008) and Wai et al. (2008). The total amount of mtDNA turnover (presented as $\sigma ={{{\displaystyle \sum }}^{\text{}}}_{i=3}^6{\nu }_i\tau {\prime }_i$, the product of turnover rate and the time for which this rate applies, summed over quiescent dynamic phases; for example, a turnover rate of 0.1 hr⁻¹ for 30 days yields σ = 0.1 × 24 × 30 = 72) is constrained more than the specific trajectory of mtDNA turnover rates, showing that a variety of time behaviours of turnover are capable of producing the observed heteroplasmy behaviour.

The bottleneck mechanism identified through our analysis has several characteristic features which facilitate experimental verification. Key among these are the prediction that heteroplasmy variance acquires an intermediate (nonzero, but not maximal) value as a result of the copy number bottleneck, then continues to increase due to mtDNA turnover in later development. Our theory also produces quantitative predictions regarding the structure of heteroplasmy distributions at arbitrary times.

The existing data that we used to perform inference and model selection display a degree of internal heterogeneity, coming from several different experimental groups. Furthermore, these data represent statistics resulting from a single pairing of mtDNA types, and it is thus arguable how conclusions drawn from them may represent the more genetically diverse reality of biology. Burgstaller et al. (2014) recently addressed this issue of a limited number of mtDNA pairings by producing novel mouse models involving mixtures of standard and several new, unexplored, wild-derived haplotypes which capture a range of genetic diversity. To test the applicability and generality of our predictions, we have perfomed new experimental measurements of germline heteroplasmy variance in these model animals under a consistent experimental protocol (see ‘Materials and methods’). We use the ‘HB’ mouse line from Burgstaller et al. (2014) pairing a wild-derived mtDNA haplotype (labelled ‘HB’ after its source in Hohenberg, Germany) with C57BL/6N; we refer to this model as ‘HB’.

Heteroplasmy measurements were taken in oocytes sampled from mice at ages 24–61 dpc (see ‘Materials and methods’ and Appendix 1; raw data in Figure 4—source data 1). The statistics of these measurements yielded $\mathbb{E}\left(h\right)$, $\mathbb{V}\left(h\right)$ and $\mathbb{V}\prime \left(h\right)$ as previously. This age range was chosen to address the regions with most power to discriminate between the competing models; the existing $\mathbb{V}\prime \left(h\right)$ data is most heterogeneous around 20–30 dpc and the later datapoints allow us to detect developmental heteroplasmy behaviour after the copy number minimum. Figure 4A shows these $\mathbb{V}\prime \left(h\right)$ measurements. The qualitative behaviour predicted by the BDP mechanism is clearly visible: variance around birth (after the copy number bottleneck) is low but non-zero, subsequently increasing with time. The ability of the BDP model to account for the magnitudes and time behaviour of heteroplasmy variance more satisfactorily than the alternative models is shown by the model fits in Figure 4A. We explored these new data quantitatively through the same model selection approach used for the existing data. As shown in Figure 4B, the BDP mechanism again experiences by far the strongest statistical support in this genetically different system.

Figure 4

Download asset Open asset

Figure 4—source data 1

Individual heteroplasmy measurements in the HB mouse model contributing to the new heteroplasmy variance data used to test our theory.

https://doi.org/10.7554/eLife.07464.007

Download elife-07464-fig4-data1-v3.zip

Figure 4C shows the result of our parameteric inference approach using these $\mathbb{V}\prime \left(h\right)$ measurements coupled with the $\mathbb{E}\left(m\right)$ measurements used previously (employing our assumption that modulation of copy number by heteroplasmy in this non-pathological haplotype is small). Strikingly, the quantitative behaviour of $\mathbb{V}\prime \left(h\right)$ with time inferred from the HB model (red) matches the previous behaviour inferred from the NZB/BALB/c system (blue) very well, suggesting that our theory is applicable across a range of genetically distinct pairings. We note that the shaded region in Figure 4C corresponds to credibility intervals around the mean behaviour of $\mathbb{V}\prime \left(h\right)$, and the fact that individual $\mathbb{V}\prime \left(h\right)$ datapoints (subject to fluctuations and sampling effects) do not all lie within these intervals is not a signal of poor model choice. An analogous situation is the observation of a scatter of datapoints outside the range of the standard error on the mean (s.e.m.), which does not imply a mistake in the s.e.m. estimate. The difference between the trace in Figure 4A and the mean curve in Figure 4C arises because Figure 4A shows the behaviour of the model under a single, optimised parameterisation, whereas Figure 4C shows the distribution of model behaviours over the posterior distributions on parameters: the mean $\mathbb{V}\prime \left(h\right)$ trace of this distribution is comparable but not equivalent to that from the single best-fit parameterisation.

To confirm more detailed predictions of our model, we also examined the specific distributions of heteroplasmy in our new measurements. Given a mean heteroplasmy and an organismal age, the parameterised BDP model predicts the structure of the heteroplasmy distribution (see ‘Materials and methods’ and next section). We parameterised the model using $\mathbb{V}\prime \left(h\right)$ values from a subset of half of the new measurements (chosen by omitting every other sampled set when ordered by time). Figure 4D shows a comparison of measured heteroplasmy distributions with a 95% bound from the parameterised BDP model. We then tested the predictions of the parameterised model against the other half of new measurements. 8 of the test measurements (2.4%) fell outside the inferred 95% bound from the training dataset, illustrating a good agreement with distributional predictions. The Anderson-Darling test was used to compare the distribution of heteroplasmy in sampled oocytes with distributions predicted by our theory (given age and mean heteroplasmy); no set of samples showed significant (p < 0.05) departures from the hypothesis that the two distributions were identical. Some example distributions are presented in Figure 4D (i), (ii), (iii).

Importantly, the BDP model yields analytic solutions for the values of all genetic properties of interest, using tools from stochastic processes (detail in ‘Materials and methods’ and Appendix 1). These results facilitate straightforward further study and fast predictions of timescales and probabilities of interest. The full theoretical approach is detailed in Appendix 1, and equations for the mean and variance of mtDNA populations and heteroplasmy are given in the Methods. In Figure 2A we illustrate that these analytic results provide an excellent match to the numeric results of stochastic simulation, a result that holds across all BDP model parameterisations. It is also straightforward to calculate the fixation probability $\mathbb{P}\left(m=0\right)$, which allows us to characterise all heteroplasmy distributions that arise from the bottlenecking process, even when highly skewed (see ‘Materials and methods’ and Appendix 1). We have thus obtained analytic solutions for the time behaviour of mtDNA copy number and heteroplasmy throughout the bottleneck with no assumptions of continuous population densities or fixed population size, under a physical model with the most statistical support given experimental data.

We can use our theory to explore the dependence of bottleneck dynamics on specific biological parameters. We first explore the effects of modulating mtDNA turnover by varying λ and ν in concert, corresponding to an increase in mtDNA degradation balanced by a corresponding increase in mtDNA replication. This increased mtDNA turnover increases the heteroplasmy variance (see Figure 5A) due to the increased variability in mtDNA copy number from the underlying random processes occurring at increased rates. We find that increasing mtDNA degradation ν without increasing λ also increases heteroplasmy variance, in addition to decreasing the overall mtDNA copy number (Figure 5B). Applying this unbalanced increase in mtDNA degradation without a matching change in replication has a strong effect on mtDNA dynamics as it corresponds to a universal change in the ‘control’ applied to the system, analogous, for example, to changing target copy numbers in manifestations of relaxed replication (Chinnery and Samuels, 1999). The simple model we use does not include feedback, and controls mtDNA dynamics solely through kinetic parameters. Perturbing the balance of these parameters thus strongly affects the expected behaviour of the system. As we discuss later, elucidation of the specific mechanisms by which control is manifest in mtDNA populations will require further research, but these numerical experiments attempt to represent the cases where a perturbation is naturally compensated for (matched changes, Figure 5A) and where it is not (unbalanced change, Figure 5B).

Figure 5

Download asset Open asset

These results suggest that an artificial intervention increasing mitochondrial degradation may generally be expected to increase heteroplasmy variance during development. An increase in mtDNA degradation is expected to either directly increase heteroplasmy variance (Figure 5B) if mtDNA populations are weakly controlled, or to provoke a compensatory, population-maintaining increase in mtDNA replication, thus increasing mtDNA turnover, which also acts to increase variance (Figure 5C) if mtDNA populations are subject to feedback control. The increase in variance through either of these pathways will increase the power of cell-level selection to remove cells with high heteroplasmy and thus purify the population. For this reason, we speculate that mitochondrial degradation may represent a potential clinical target to address the inheritance of mtDNA disease (more detail in Appendix 1).

Our model also allows us to explore the effect of different mtDNA types experiencing different selective pressures, by setting λ₁ ≠ λ₂ (mutant mtDNA experiences a proliferative advantage or disadvantage). Such a selective difference causes changes in both mean heteroplasmy and heteroplasmy variance, as shown in Figure 5C (e.g., if heteroplasmy decreases towards zero, heteroplasmy variance will also decrease, as the wildtype is increasingly likely to become fixed). We do not focus further on selection in this study, noting that selective pressures are likely to be specific to a given pair or set of mtDNA types and are not generally characterised well enough to perform satisfactory inference. However, we do note that our theory gives a straightforward prediction for the functional form of mean heteroplasmy when nonzero selection is present, a sigmoid with slope set by the fitness difference (see ‘Materials and methods’).

A key feature of mtDNA diseases is that pathological symptoms usually manifest when heteroplasmy in a tissue exceeds a certain threshold value, with few or no symptoms manifested below this threshold (Rossignol et al., 2003). The probability and timescale with which cellular heteroplasmy may be expected to exceed a given value is thus a quantity of key interest in clinical planning of mtDNA disease strategies.

In our model, the probability, as a function of time, of a cell containing m₁ wildtype and m₂ mutant mtDNAs can be straightforwardly derived. The resultant analytic expression involves a hypergeometric function, also an important mathematical element in expressions describing mtDNA statistics based on classical population genetics (Kimura, 1955; Wonnapinij et al., 2008). The probability of obtaining a given heteroplasmy can therefore be computed as a sum over all copy number states that correspond to that heteroplasmy. However, as hypergeometric functions are comparatively unintuitive and computationally expensive, we here employ an approximation to the distribution of heteroplasmy based upon the above moments that are straightforwardly calculable from our model. This approximation involves fixation probabilities for each mtDNA type and a truncated Normal distribution for intermediate heteroplasmies (see ‘Materials and methods’). In Appendix 1 we show that this approximation corresponds well to the exact distributions calculated using the hypergeometric function. We underline that exact heteroplasmy distributions are straightfoward to compute using our approach: we use the truncated Normal approximation as it represents the exact distribution well, is more intuitively interpretable, and is computationally very inexpensive.

Using this approach, the probability with time of a cell exceeding a threshold heteroplasmy h^* can be straightforwardly computed for any initial heteroplasmy, allowing rigorous quantitative elucidation of this important clinical quantity (see ‘Materials and methods’). Figure 5D illustrates this computation by showing the analytic probability with which thresholds h^* = 0.4, 0.5, 0.6 are exceeded at a time t, given the example initial heteroplasmy h = 0.3. These results serve as a simple example of the power of our modelling approach: any other specific case can readily be addressed. Our theory thus allows general quantitative calculation of the probability (and timescale) that any given heteroplasmy threshold will be exceeded, given knowledge of the initial (or early) heteroplasmy.

We next turn to the question of estimating heteroplasmy levels in a developed organism by sampling cells during development. This principle, clinically termed preimplantation genetic diagnosis (Steffann et al., 2006; Poulton et al., 2009), assists in clinical planning by allowing inference of the specific heteroplasmic nature of the embryo itself rather than a population average of an affected mother's oocytes (Treff et al., 2012). However, the complicated and stochastic nature of the bottleneck makes this inference a challenging problem.

Given a heteroplasmy measurement from sampling h_m, accurate preimplantation diagnosis is contingent on knowledge of the distribution $\mathbb{P}\left(h|h_m\right)$, that is, the probability that the embryonic heteroplasmy is h given that a measurement h_m has been made. We can use our modelling framework and Bayes' theorem (see ‘Materials and methods’) to obtain a formula for this conditional probability, allowing a rigorous probability to be assigned to inferences from preimplantation sampling. Here, as above, we employ the truncated Normal approximation for the heteroplasmy distribution, noting that the exact treatment using hypergeometric functions is straightforward but more computationally expensive. Figure 5E illustrates this process by showing the probability distributions on embryonic heteroplasmy when measurements h_m = 0.1 or 0.4 have been taken at different times during development. The increasing heteroplasmy variance through development means that substantially greater uncertainty is associated with heteroplasmy values inferred using measurements taken at later times. In conclusion, although care must be taken in applying this reasoning to cell types in which, for example, mitochondrial and cell turnover rates differ from those assumed here, or differentiation leads to tissue-specific selective factors acting on the mtDNA population, this formalism provides a general means of rigorously inferring embryonic heteroplasmy through genetic diagnosis sampling.

Our ‘bottom-up’ model represents individual mtDNAs as elements which replicate and degrade either randomly or deterministically according to the model parameterisation. Consonant with experimental studies showing that it is often a single mutant genotype that dominates the non-wildtype mtDNA population of a cell (Khrapko et al., 1999), we consider two mtDNA types (wildtype and mutant), though our model can readily be extended to more mtDNA types. We denote the number of ‘wild-type’ mtDNAs in a cell as m₁ and the number of ‘mutant’ mtDNAs as m₂. The heteroplasmy of a cell is then $h=\frac{m_2}{m_1+m_2}$, that is, the population proportion of mutant mtDNA.

Individual mtDNAs are capable of replication and degradation, with rates denoted λ and ν respectively. According to a binary categorical parameter S, these events may be deterministic (S = 0; the mtDNA population replicates and degrades by a fixed amount per unit time) or Poisson processes (S = 1; each individual mtDNA randomly replicates and degrades with average rates λ and ν). A parameter α controls the proportion of mtDNAs capable of replication: α = 1 allows all mtDNAs to replicate throughout development, α < 1 enforces a subset proportion α of replicating mtDNAs a time cutoff T after conception.

A parameter c (cluster size; a non-negative integer) dictates the partitioning of mtDNAs at cell divisions. When c = 0, partitioning is deterministic, so each daughter cell receives exactly half of its parent's mtDNA. For c > 0, partitioning is stochastic. When c = 1, partitioning is binomial: each mtDNA has a 50% chance of being inherited by either daughter cell. When c > 1, the parent cell's mtDNAs are grouped in clusters of size c before division. Each cluster is then partitioned binomially, with a 50% chance of being inherited by either daughter cell.

The mtDNA population changes in different ways as development progresses, first decreasing, then recovering, then slowly growing. We include the possibility of different ‘phases’ of mtDNA dynamics in our model to capture this behaviour. Each phase j has its own associated pairs of λ_j, ν_j parameters and may either be quiescent (involving no cell divisions) or cycling (encompassing n_j cell divisions). Thus, we may have an initial cycling phase with low mtDNA replication rates, so that copy number falls for several cell divisions, then a subsequent ‘recovery’ cycling phase with higher replication rates so that mtDNA levels are amplified, then quiescent phases as cell lineages are identified. We allow six different phases, with the first two fixed as cycling phases with the above doubling times, and the final phase fixed to include no mtDNA replication (representing the stable, final occyte state).

The initial conditions of our model involve an initial mtDNA copy number m₀ (the total number of mtDNAs in the fertilised oocyte) and an initial heteroplasmy h₀ (the fraction of these mtDNAs that are mutated).

We used three datasets for mtDNA copy number during mouse development: Cao et al. (2007); Cree et al. (2008); and Wai et al. (2008). We use two datasets for heteroplasmy variance during development: Wai et al. (2008) and Jenuth et al. (1996). By convention, we use the normalised versions of heteroplasmy variance (i.e., measured variance divided by a factor h(1 − h)). Where the measurements were not given explicitly in these publications, we manually analysed the appropriate figures to extract the numerical data. For these values, we used data from correspondence regarding the Wai study (reply to [Samuels et al., 2010]), and manually normalise the Jenuth dataset. The Jenuth dataset contains measurements taken in ‘mature oocytes’ with no time given; we assume a time of 100 dpc for these measurements, though this time is generalisable and does not qualitatively affect our results. All values are presented in Appendix 1. Data on cell doubling times in the mouse germ line is taken from Lawson and Hage (1994), suggesting that doubling times start with an interval of every 7 hr, then after around 8.5 days post conception (dpc) increase to 16 hr, before the onset of a quiescent regime around 13.5 dpc (roughly consistent with the estimate of ∼25 divisions between generations in the female mouse germ line [Drost and Lee, 1995]).

We use Gillespie algorithms, also known as stochastic simulation algorithms (Gillespie, 1977), to explore the behaviour of our model of the bottlenecking process for a given parameterisation. For a given model parameterisation, the Gillespie algorithm is used to simulate an ensemble of 10³ possible realisations of the time evolution of mtDNA content, and the statistics of this ensemble are recorded. The experimental data we use is derived from sets of measurements of different sizes; to compare simulation data with an experimental datapoint i corresponding to a statistic derived from n_i measurements, we sampled a random subset of n_i of the 10³ simulated trajectories (all datapoints but one have n ≪ 10³), and used this subset to derive the simulated statistic for comparison to datapoint i (Johnston, 2014).

To fit the different models to experimental data we define a distance measure, a sum-of-squares residual between the $\mathbb{E}\left(m\right)$ (in log space) and $\mathbb{V}\left(h\right)$ dynamics produced by our model and observed in the data, weighted to facilitate comparison of these different quantities (Johnston, 2014). We also constrain copy number to be <5 × 10⁵ at all points throughout development, rejecting parameterisation that disobey this criterion. Metropolis MCMC was used to identify the best-fit parameterisation according to this distance function. For statistical inference, we use approximate Bayesian computation (ABC), a statistical approach that has successfully been applied to parametric inference and model selection in dynamical systems (Toni et al., 2009) to infer posterior probability distributions both for individual models and the parameters of the models given experimental data. ABC samples posterior probability distributions on parameters that lead to behaviour within a certain threshold distance of the given data; these posteriors are shown to converge on the true posteriors as the threshold value decreases to zero (see Appendix 1). We employed an MCMC sampler with randomly-selected initial conditions. For further details, including priors, thresholds and step sizes used in ABC, see Appendix 1. Minimum copy number was recorded directly from the resulting trajectories; our measure of total turnover σ is defined as $\sigma ={{{\displaystyle \sum }}^{\text{}}}_{i=3}^6\tau {\prime }_i{\nu }_i$, the sum over quiescent dynamic phases of the product of degradation rate and phase length.

Heteroplasmic mice were obtained from a heteroplasmic mouse line (HB) we created previously by ooplasmic transfer (Burgstaller et al., 2014). This mouse line contains the nuclear DNA of the C57BL/6N mouse, and mtDNAs both of C57BL/6N and a wild-derived house mouse. Both mtDNA variants belong to the same subspecies, Mus musculus domesticus. For details on sequence divergence (see Burgstaller et al., 2014).

Mice were sacrificed at the indicated ages by cervical dislocation. Ovaries were extracted and immediately placed in cryo-buffer containing 50% PBS, 25% ethylene glycol and 25% DMSO (Sigma–Aldrich, Austria) and stored at −80°C. For oocyte extraction, ovaries were placed into a drop of cryo-buffer and disrupted using scalpel and forceps. Oocytes were collected and remaining cumulus cells were removed mechanically by repeated careful suction through glass capillaries. Prepared oocytes were then washed in PBS before they were individually placed into compartments of 96-well PCR plates (Life Technologies, Austria) containing 10 μl of oocyte-lysis buffer (Lee et al., 2012) (50 mM Tris-HCl, [p.H 8.5], 1 mM EDTA, 0.5% tween-20 [Sigma–Aldrich, Austria] and 200 μg/ml Proteinase K [Macherey–Nagel, Germany]). Samples covered stages from primary oocytes of 3 day-old mice up to mature oocytes of 40 day-old mice. Samples were lysed at 55°C for 2 hr, and incubated at 95°C for 10 min to inactivate Proteinase K. The cellular DNA extract was finally diluted in 190 μl Tris-EDTA buffer, pH 8.0 (Sigma–Aldrich, Austria). 3 μl were used per qPCR reaction.

Heteroplasmy quantification was performed by ARMS-qPCR, an established method in the field (Steinborn et al., 2000; Paull et al., 2013; Tachibana et al., 2013), as described in Burgstaller et al. (2014). The study was conducted according to MIQE (minimum information for publication of quantitative real-time PCR experiments) guidelines (Bustin et al., 2009; Burgstaller et al., 2014). The proportion between HB derived and C57BL/6N mtDNA was determined by ARMS-qPCR assays based on a SNP in mt-rnr2 (Burgstaller et al., 2014). These assays were normalised to changes in the input mtDNA amount by consensus assays, located in conserved regions of mt–Co2 and mt–Co3 (see Appendix 1). For the calculation of mtDNA heteroplasmy, the assay detecting the minor allele (C57BL/6N or wild-derived <50%) was always used. If both specific assays gave values >50% (which can happen around 50% heteroplasmy), the mean value of both assays was taken. All qPCR runs contained no template controls (NTCs) for all assays; these were negative in 100%. Further experimental details available in Appendix 1.

In the BDP model, processes within a cell cycle constitute a birth-death process which can be solved using generating functions (Gardiner, 1985). For binomial partitioning, the generating function for the system after an arbitrary number of divisions has a recursive structure (Rausenberger and Kollmann, 2008; Johnston and Jones, 2015) and an analytic solution can be obtained through solving a Riccati recurrence relation. This reasoning also extends to the different phases of replication and degradation, allowing an exact generating function to be constructed for an arbitrary point in the bottleneck. Derivatives of this generating function are then used to obtain moments of the distributions of interest. The full procedure is given in Appendix 1. Recall that we assume that the bottlenecking process consists of a series of dynamic phases, which may either involve cycling cells (and hence cell divisions) or quiescent cells. The expression for mean mtDNA copy number $\mathbb{E}\left(m,t\right)$ at time t is:

where n_i is the number of cell divisions in phase i (0 for quiescent phases), τ_i is the length of a cell cycle in cycling phase i, $\tau {\prime }_i$ is the time spent in quiescent phase i (0 for cycling phases), and ${\tau }^{\text{*}}={\text{Σ}}_i\left(n_i{\tau }_i+\tau {\prime }_i\right)$, so that t − τ^* is the time since the last cell division. $\mathbb{E}\left(m,t\right)$ is thus intuitively interpretable as a product of the initial copy number with the effects of halving at each cell division, and the copy number evolution through past and current cell cycles and quiescent phases.

The expression for the variance is lengthier, taking the form

where Φ is a lengthy, though algebraically simple, function of all physical parameters, which we derive and present in Appendix 1. Once the means and variances associated with mutant and wild-type mtDNAs have been determined (for brevity, we write these as ${\mu }_1\equiv \mathbb{E}\left(m_1,t\right),{\sigma }_1^2\equiv \mathbb{V}\left(m_1,t\right)$ and ${\mu }_2\equiv \mathbb{E}\left(m_2,t\right),{\sigma }_2^2\equiv \mathbb{V}\left(m_2,t\right)$), the relations below can be used to compute heteroplasmy statistics:

The predicted mean heteroplasmy at time t assuming a constant selective pressure (though this assumption can straightforwardly be relaxed) is given by Equation 4, which, given Equation 2, straightforwardly reduces to

where h₀ is initial heteroplasmy and Δλ is the increase (or decrease, if negative) in replication rate of mutant over wild-type mtDNA. Equation 6 predicts that mean heteroplasmy in the presence of selection will follow a sigmoidal form (as expected from population dynamics [Futuyma, 1997], by the constraint that h₀ must lie between 0 and 1, and by the intuitive fact that heteroplasmy changes slow down as these limits are approached).

The probability of heteroplasmy exceeding a certain threshold h^* is simply given by integrating the probability distribution of heteroplasmy between h^* and 1. The exact distribution of heteroplasmy can be written as a sum over hypergeometric functions; however, for computational efficiency and interpretability, we employ an approximation to this distribution involving the truncated Normal distribution and fixation probabilities. As shown in Appendix 1, the distribution of heteroplasmy, taking possible fixation into account, can be well approximated by

where $\mathcal{N}\prime$ is the truncated Normal distribution (truncated at 0 and 1), μ and σ² are found numerically given our model results for $\mathbb{E}\left(h\right)$ and $\mathbb{V}\left(h\right)$, and ${\zeta }_1\equiv \mathbb{P}\left(h=0\right)$ and ${\zeta }_2\equiv \mathbb{P}\left(h=1\right)$ are fixation probabilities, also straightforwardly calculable from our model. The probability of threshold crossing for 0 < h^* < 1 is then

Given a sampled measurement heteroplasmy h_m, the probability $\mathbb{P}\left(h_0|h_m\right)$ that embryonic heteroplasmy is h₀ is given by Bayes' theorem $\mathbb{P}\left(h_0|h_m\right)=\mathbb{P}\left(h_m|h_0\right)\mathbb{P}\left(h_0\right)/\mathbb{P}\left(h_m\right)$. Assuming a uniform prior distribution on embryonic heteroplasmy (though this can be straightforwardly generalised), we thus obtain $\mathbb{P}\left(h_0|h_m\right)=\mathbb{P}\left(h_m|h_0\right)/{{{\displaystyle \int }}^{\text{}}}_0^1\mathbb{P}\left(h_m|h_0\prime \right)d{h}_0\prime$, and using the above expression for the heteroplasmy,

where μ, σ², ζ₁, ζ₂ are functions of h₀: μ, σ² may be found numerically and the ζ values are analytically calculable (see Appendix 1).

Table 1 contains the datapoints used in this study. These data are taken from Tables 1, 2 of Cree et al. (2008) (labelled ‘Cao’); Tables 1, 2 of Cao et al. (2007); Appendix figure 1 of Wai et al. (2008) and Appendix figure 1 of the following correspondence (reply to Samuels et al., 2010) (labelled ‘Wai’); and Table 2 of Jenuth et al. (1996) (labelled ‘Jenuth’). Convention in the literature suggests that normalisation of measured heteroplasmy variance values, performed by division by a factor h(1 − h), allows comparison of variance values from lines with diverse absolute heteroplasmies: the Wai data from correspondence is already normalised, and we manually normalised the Jenuth data using the h values present.

Appendix figure 1

Download asset Open asset

Where data in the original studies were presented as a function of number of cells in a developing organism, as opposed to an explicit function of time, we have assigned times using the 7 hr → 16 hr doubling times from Lawson and Hage (1994). Other sources assume a 15 hr doubling time throughout early development: using the data interpreted in this way did not lead to a qualitative difference in our conclusions and very little quantitative change in posterior distributions (data not shown). Some datapoints did not have associated or readily available sample sizes N: for these datapoints we estimated N using available evidence in the publication. To check for dependence on these values of N we performed our inference process with a range of alternative N values and with a test case where N was set to 100 for every datapoint: all results and posteriors were qualitatively similar, showing a lack of strong dependence of our conclusions on the specific numbers of samples involved in deriving the experimental measurements (data not shown).

The specific units of inheritance of mtDNA have been debated in the literature for decades. The smallest possible unit of inheritance is a single mtDNA molecule; some studies have hypothesised that the unit of inheritance consists of groups of mtDNA molecules. Within this picture, debate exists as to whether these groups are semi-permanent associations of molecules (which we will refer to as ‘quenched’ sets) or more fluid transient colocalisations of molecules (which we refer to as ‘unquenched’). Furthermore, the size of these units is debated, with estimates ranging from an average size of 1.4–10 mtDNA molecules (Bogenhagen, 2012; Kukat and Larsson, 2013), and it is unknown whether the mtDNAs within a group are strictly homoplasmic or if heteroplasmic groups are possible, although current evidence, at the finest resolution, points towards homoplasmic groups of size <2 (Gilkerson et al., 2008; Poe et al., 2010; Wallace and Chalkia, 2013; Jakobs and Wurm, 2014).

We will classify these different pictures with three parameters. First, the characteristic size c of an mtDNA group. Second, a classifier denoting whether these groups are quenched (in the sense that the individual constituents of a group remain the same over many cell divisions) or unquenched (in the sense that the individual constituents of a group may change between cell divisions). Third, a classifier denoting whether groups are necessarily homoplasmic, or if heteroplasmy is permitted.

An early hypothesis from Jacobs et al. (2000) considered ‘nucleoids’ which correspond to quenched heteroplasmic groups with c > 1, retaining their internal structure across cell divisions and containing different mtDNA types. If the mitochondrial organelle is the unit of inheritance, we may expect unquenched heteroplasmic groups with c > 1, as mitochondrial dynamics act to mix the content of the mitochondrial system between cell divisions, but organelles are likely to contain more than one mtDNA molecule. If nucleoids are the units of inheritance and, as current understanding suggests, nucleoids are small and homoplasmic (if mtDNA indeed exists in groups at all), the appropriate picture is c ≃ 1, homoplasmic groups.

Here we show that the heteroplasmy statistics resulting from these different pictures of grouped inheritance collapse onto two representative cases: first, that corresponding to homoplasmic clusters with c > 1, and second, that corresponding to c = 1 (binomial inheritance). Quenching—whether mtDNA content can remix within nucleoids—is shown to be unimportant in determining heteroplasmy statistics. Our model for these different situations is as follows. We consider a cellular population as consisting of a set of mtDNA molecules, existing in groups of size c. During a cell cycle, the population of groups doubles deterministically (we ignore random birth-death dynamics in this model, in order to focus on partitioning dynamics), so that every group produces one exact copy of itself. For unquenched simulations, a new set of groups is then formed by resampling the individual mtDNA constituents of the cell. For quenched simulations (representing the situation postulated in Jacobs et al. (2000)), the existing groups remain intact. At cell divisions, groups are binomially partitioned between the two daughter cells.

The model is initialised with a cell containing m₀ mtDNAs, split into (1 − h)m₀ wild-type and hm₀ mutant molecules. These mtDNAs are clustered into m₀/c groups, according to the cluster picture under consideration (i.e., homoplasmic or heteroplasmic clusters). We simulate the subsequent doubling then partitioning of this system through cell divisions many times, assuming a constant cell cycle length, and record the cell-to-cell heteroplasmy variance with time.

Appendix figure 1 shows the resultant heteroplasmy variance trajectories for different cases (with h₀ = 0.1; other initial heteroplasmies showed similar behaviour). The first striking result is that the inheritance of heteroplasmic groups produces the same heteroplasmy variance as binomial partitioning, regardless of cluster size. This behaviour is due to the balance between stochasticity associated with the makeup of, and partitioning of, groups. A small number of large groups will experience substantial partitioning noise, but larger heteroplasmic groups are more likely to faithfully represent the overall cell heteroplasmy. As identified in Jacobs et al. (2000), the inheritance of heteroplasmic groups thus provides a means to facilitate local mtDNA complementation while provoking no increase in heteroplasmy variance beyond that associated with binomial partitioning of elements at divisions.

We also observe that quenched populations behave in the same way as unquenched populations. In the case of homoplasmic groups, this result is obvious, as a set of homoplasmic nucleoids of a given size can only be constructed in one way for a given number of mtDNA molecules of different types. For heteroplasmic groups, this result implies that resampling the cellular population to produce a new group produces a negligible amount of additional stochasticity compared to that already present in the random makeup and inheritance of groups. Thus, the only determinant factors of heteroplasmy variance related to the inheritance of groups are whether groups are homoplasmic or heteroplasmic, and, if the former, the characteristic size of groups.

These results illustrate that the binomial inheritance model can also describe the statistics associated with heteroplasmic nucleoids of arbitrary size, over a timescale of several dozen cell divisions (suitable to describe the developmental process). The theoretical long-term behaviour of these systems involves some more subtleties. At much longer times, the probability that all mtDNA types become extinct in a cell is not negligible. When complete extinction cannot be ignored, heteroplasmy statistics become poorly defined. This extinction timescale is shorter for cluster inheritance than for binomial inheritance, as a greater variability in copy number (though not in heteroplasmy) results from each division for larger clusters. However, our simulations indicate that as long as the heteroplasmy variance associated with heteroplasmic clusters remains well defined, it matches that resulting from binomial inheritance.

We propose that a reasonable view may be that individual mitochondrial fragments, including several small, homoplasmic nucleoids, are the likely elements of inheritance at partitioning. Furthermore, there is likely some movement of these nucleoids within the mitochondrial network, and fission and fusion likely mean that a given mtDNA will not be associated with the same static mitochondrial element in perpetuity. In this case, the picture of an unquenched, heteroplasmic group of mtDNAs—those contained within a discrete element of the mitochondrial system—seems most reasonable. We can thus speculate that, as demonstrated by the previous results, the precise size of mitochondrial fragments at partitioning is not important for the heteroplasmy dynamics (nor indeed is whether they are quenched or unquenched). Our simple binomial partitioning model is thus consistent with what one might consider the most physiologically plausible model, and indeed with any models not involving large and strictly homoplasmic groups as the elements of mitochondrial inheritance.

Our model is a function of the parameter set θ = {ν_i, λ_i, n_i, τ_i, S, α, T, c, h₀, m₀, δλ}. For reference, the meanings of these parameters are (as in Figure 1D in the Main text): replication (λ_i) and degradation (ν_i) rates; number of cell divisions (n_i) and cell cycle length (τ_i) in each dynamic phase i; deterministic or stochastic dynamics label (S = 0, 1 respectively); a proportion α of mtDNAs capable of replication after threshold time T; deterministic (c = 0), binomial (c = 1) or clustered (c > 1) partitioning at divisions; initial heteroplasmy h₀ and initial copy number m₀; δλ is an additional parameter allowing a possible difference in replicative rates between mutant and wildtype mtDNA: this is zero unless otherwise stated. For the following parameters we use uninformative uniform priors on the given interval: λ_i, ν_i ∈ [0,1] hr⁻¹; S ∈ {0,1}; α ∈ [0.005,1]; T ∈ [0,100] day; c ∈ [0,100] day; h₀ ∈ [0,1]; m₀ ∈ [0,10⁶]. The following values are fixed from experimental studies (Lawson and Hage, 1994; Drost and Lee, 1995): n₁ = 29; n₂ = 7; τ₁ = 7 hr; τ₂ = 16 hr. λ₆ = 0 hr⁻¹ is fixed to avoid mtDNA proliferation after development; h₀ = 0.2 is fixed as an intermediate value as heteroplasmy variance measurements are generally normalised; δλ, a parameter allowing a difference in replicative rates between mutant and wildtype mtDNAs, is fixed at zero throughout as we ignore selective pressure. The parameter τ_i for i > 2 is used to determine the length of time spent in different quiescent phases and is subject to the uniform prior τ_i ∈ [0, 50] day.

Given these priors, we use an approximate Bayesian computation (ABC) approach to build a posterior distribution over the parameters in our bottlenecking model (Toni et al., 2009). ABC involves using a summary statistic $\rho \left(\theta ,\mathcal{D}\right)$ to compare the available data $\mathcal{D}$ to the predictions of a model given parameters θ. If parameter sets are sampled from the set for which ρ ≤ ϵ, where ϵ is a threshold difference between the resulting model behaviour and experimental data, the posterior distribution $P\left(\theta |\rho \left(\theta ,\mathcal{D}\le \mathit{ϵ}\right)\right)$ is sampled, which is argued to sufficiently approximate $P\left(\theta |\mathcal{D}\right)$ for suitably small ϵ (Marjoram et al., 2003).

We define a residual sum-of-squares difference between the results of a simulated model and experimental data (Johnston, 2014):

where $\mathcal{D}$ denotes experimental data. We thus amalgamate experimental results of two types: mean mtDNA copy number (with N_m data points measuring ${\mathbb{E}}_{\mathcal{D}}\left(m\right)$ at times $t_{\mathcal{D},m}^{\left(i\right)}$), and mean and variance of heteroplasmy (with N_h data points measuring ${\mathbb{V}}_{\mathcal{D}}\left(h\right)$ at times $t_{\mathcal{D},h}^{\left(i\right)}$). The sets of data for $\mathbb{E}\left(m\right)$ and $\mathbb{V}\left(h\right)$ contain different numbers of points and are of different absolute magnitudes. We compensate for these differences by using the logarithms of copy number measurements (as these values span several orders of magnitude), and weighting parameter A₁ = 10³. This weighting parameter compensates for the different magnitudes and number of datapoints in each class of measurement, ensuring that the contribution to the total residual from each set of data is of comparable magnitude. Our summary statistic thus records a residual sum-of-squares difference between experiment and simulation values for $\text{log}\mathbb{E}\left(m\right)$ and $\mathbb{V}\left(h\right)$ at each time point where an experimental measurement exists.

We performed our model selection process using several different alternative protocols, including comparing logarithms of $\mathbb{V}\left(h\right)$ measurements (in contrast to the raw values) and varying A₁ over orders of magnitude from 10²–10⁴ (corresponding to unbalanced weighting, favouring $\mathbb{E}\left(m\right)$ and $\mathbb{V}\left(h\right)$ data respectively). In all cases, the BDP model identified in the Main text experienced substantially more support than any alternative. For inference involving the new dataset from the HB model system, we use the default protocol above and set A₁ = 3 × 10³ to account for the threefold decrease in available $\mathbb{V}\prime \left(h\right)$ datapoints.

We use an MCMC implementation of ABC, whereby we construct a Markov chain θ_i, where each state consists of a set of trial parameters to be assessed. We create θ_i+1 by perturbing each parameter within θ_i with a perturbation kernel consisting of a Normal distribution on each parameter with standard deviations between 0.1–1% of the width of the prior (varied as the model depends more sensitively on some parameters than others). In the case of discrete parameter c, a continuous representation c′ is used and varied in the MCMC approach, with $c=\lfloor 100c\prime \rfloor$. We accept θ_i+1 as the new state of the chain if $\rho \left({\theta }_{i+1},\mathcal{D}\right)\le \mathit{ϵ}$. We ran 10⁶ MCMC iterations for ABC model selections and checked convergence by running five instances of each simulation for different random number seeds.

For the initial optimisation of model fitting, we ran 10⁶ MCMC steps using the protocol above but accepting a move according to the Metropolis–Hastings protocol (Hastings, 1970), recording the parameterisation leading to the lowest recorded residual. In this case we used uninformative initial conditions, with identical choices for all rate parameters, corresponding to an inaccurate trajectory of copy number and heteroplasmy variance. For model selection, we used the protocol above, with a different set of parameters θ^M for each model M, with each MCMC step proposing a random model from the Cao alone, Wai alone, and BDP set described in the text, as in the SMC ABC model selection protocol proposed in Toni et al. (2009). We record the proportion of accepted steps involving each model type. The parameterisations found through initial optimisation were used as initial conditions in the ABC model selection and inference simulations.

Initial optimisation identified parameterisations all displaying residuals under ϵ = 50. We chose ϵ = {45, 50, 60, 75} for the ABC model selection simulations to display the varying degrees of support for each model as stricter agreement with experiment was enforced. We chose ϵ = 40 for the ABC inference of BDP model parameterisation to ensure these models all displayed better fits to data than the alternative models. In Appendix figure 2 we illustrate the distribution of squared residuals for the BDP model under a range of ϵ values.

Appendix figure 2

Download asset Open asset

In Appendix figure 3 we display all posterior distributions for all parameters resulting from our ABC approach assuming the BDP model. There is substantial variability in the possible timescales and magnitudes of random turnover associated with each random dynamic phase i > 2, exemplified by the complicated and bimodal structure of the posteriors on these parameters. This variability reflects the fact that an increase in heteroplasmy variance can be achieved through a variety of specific mtDNA trajectories, and current experimental data is insufficient to distinguish specific time behaviours within this variety. However, the total contribution of each random phase to the overall dynamics is more constrained, as shown in the posterior distribution on a measure of total random turnover $\sigma ={{{\displaystyle \sum }}^{\text{}}}_{i=3}^6\tau {\prime }_i{\nu }_i$. This quantity is the sum over all later phases of the product of the length of that phase and the rate of random turnover, thus giving a measure of total random turnover. The fact that this posterior is more tightly constrained than the posteriors on individual t_i, ν_i parameters suggests that the required mtDNA turnover can be achieved through a range of specific dynamic trajectories from the inferred mechanism: for example, the exact time at which random mtDNA turnover sharply increases is currently flexible (though constrained to lie around 25 dpc) without more detailed data. This flexibility is also observed in the trajectories of posterior distributions in the Main text.

Appendix figure 3

Download asset Open asset

Table 2 contains the measurements of heteroplasmy h, mean heteroplasmy $\mathbb{E}\left(h\right)$, raw and normalised heteroplasmy variance $\mathbb{V}\left(h\right)$ and $\mathbb{V}\prime \left(h\right)$, and number of datapoints n, from the HB model system. Experimental procedures are described in ‘Materials and methods’; further specifics follow.

Consensus assays:

Co2-f: GCCAATAGAACTTCCAATCCGTATAT,

Co2-r: TGGTCGGTTTGATGTTACTGTTG,

Co2-FAM: CTGATGCCATCCCAGGCCGACTAA-BHQ1 (Amplicon length: 136 bp);

Co3-f: TCTTATATGGCCTACCCATTCCAA,

Co3-r: GGAAAACAATTATTAGTGTGTGATCATG,

Co3-FAM: TTGGTCTACAAGACGCCACATCCCCT-BHQ1 (Amplicon length: 103 bp).

ARMS-assays:

16SrRNA2340(3)G-f: AATCAACATATCTTATTGACCaAG (haplotype C57BL/6N),

16SrRNA2340(3)A-f: AATCAACATATCTTATTGACCgAA (haplotype HB);

16SrRNA2458-r: CAC CAT TGG GAT GTC CTG ATC,

16SrRNA-FAM: FAM-CAA TTA GGG TTT ACG ACC TCG ATG TT-BHQ-1. (Amplicon length: 142 bp).

Every qPCR run consisted of one consensus and an ARMS assay.

Master-mixes for triplicate qPCR reactions contained 1× buffer B (Solis BioDyne, Estonia); 4.5 MgCl₂ for the ARMS and the Co3 consensus assays, and 3.5 mM MgCl₂ for the Co2 consensus assay; 200 M of each of the four deoxynucleotides (dNTPs, Solis BioDyne, Estonia), HOT FIREPol DNA polymerase according to the manufacturers instructions (Solis BioDyne, Estonia), 300 nM of each primer and 100 nM hydroloysis probe (Sigma–Aldrich, Austria). Per reaction 12 μl of master-mix and 3 μl DNA were transferred in triplicates to 384-well PCR plates (Life Technologies, Austria) using the automated pipetting system epMotion 5075TMX (Eppendorf, Germany). Amplification was performed on the ViiA 7 Real-Time PCR System using the ViiA 7 Software v1.1 (Life Technologies, USA). DNA denaturation and enzyme activation were performed for 15 min at 95°C. DNA was amplified over 40 cycles consisting of 95°C for 20 s, 58°C for 20 s and 72°C for 40 s for both assays.

The standard curve method was applied. Amplification efficiencies were determined for each run separately by DNA dilution series consisting of DNA from wild-derived mice, harbouring the respective analysed mtDNA. Typical results: slope = −3.665, −3.562, −3.461, −3.576; mean efficiency = 0.87, 0.9, 0.94, 0.90; and Y-intercept = 32.2, 28.3, 33.8, 34.5; for the consensus Co2, consensus Co3, C57BL/6N and HB assays respectively. Coefficient of correlation was >0.99 in all assays in all runs. All target samples lay within the linear interval of the standard curves. To test for specificity, in each run a negative control sample, that is, a DNA sample of a mouse harbouring the mtDNA of the non-analysed type in the heteroplasmic mouse (i.e., C57BL/6N or HB mtDNA) was measured. All assays could discriminate between C57BL/6N and HB mtDNA at a minimum level of 0.2%. Target sample DNA was tested for inhibition by dilution in Tris-EDTA buffer (Sigma–Aldrich, Austria), pH 8.0.

Here we summarise potential mechanisms for the bottleneck that conflict with our statistical interpretation, highlighting the reasons for the conflict. We also propose further experiments that would efficiently provide more evidence to distinguish these hypotheses.

Wai et al. (2008) performs BrU labelling to observe the proportion of mitochondria replicating in primary oocytes between P1-4 (21–25 dpc on our time axis). The observations contained therein (Appendix figure 2 in Wai et al., 2008) show a small subset of BrU-labelled mitochondrial foci compared to the overall population of mitochondria labelled with another dye. Here we show that this observation is compatible with (and expected from) our proposed model of random mtDNA turnover.

Consider a population of mtDNAs replicating with rate λ and degrading with rate ν. We model the BrU labelling assay as follows. At time t = 0, we begin the BrU labelling, which we conservatively model as a perfect process, so that every mtDNA that replicates becomes labelled. We continue this labelling until t = t^*, when we observe the proportion of labelled mtDNAs.

For simplicity, we will consider a fixed population of mtDNAs of size N, though this reasoning extends to changing population size. We denote by l the number of labelled mtDNAs. After BrU exposure, this number may change in three ways: (A) a replication event involving a previously unlabelled mtDNA will produce two new labelled mtDNAs; (B) a replication event involving a previously labelled mtDNA will produce one new labelled mtDNA; (C) a degradation event involving a labelled mtDNA will remove one labelled mtDNA. The dynamics of labelled mtDNA number during BrU exposure are given by

Assuming that l = 0 at t = 0, the solution of this equation, for the number of labelled mtDNAs at time t^*, is

Assuming a constant population size requires λ = ν. The conclusions of this illustrative study do not substantially change if we allow (λ ≠ ν) and hence an increasing or decreasing population. We consider the values of λ and ν required to yield values of l comparable with those found in Wai et al. (2008). We will roughly estimate these values, based on the proportion of labelled foci observable, as l = 0.5 N for 24 hr BrU exposure (half of observed mtDNAs being labelled) and l = 0.05 N for 2 hr BrU exposure (5% of observed mtDNAs being labelled). A value of λ = ν = 0.014 hr⁻¹ yields l = 0.49 N at t^* = 24 hr and l = 0.055 N at t^* = 2 hr.

Figure 3 in our Main text gives the posterior distribution on ν, characterising the rate of random mtDNA turnover in our model, at different times. It can be seen that a value of ν = 0.35 day⁻¹ comfortably falls within the region of high posterior density during the time range 21–25 dpc—lying immediately before the strong increase in random turnover that our model subsequently predicts. Our inferred mechanism of random mtDNA turnover is thus compatible with the observations of a labelled subset of mtDNAs in the BrU incorporation assay in Wai et al. (2008)—we would expect to see roughly the observed labelling proportion simply due to the likely rates of random mtDNA turnover inferred at that stage of development. Furthermore, we can use this line of reasoning to produce a testable prediction: similar experiments carried out several days later—when random mtDNA turnover is inferred to increase substantially—should show a larger subset of labelled mtDNAs for the same BrU exposure.

The results from our model suggest a potential clinical pathway for increasing heteroplasmy variance, and thus the power of the bottleneck to remove heteroplasmic cells. We have shown that upregulation of mtDNA degradation (e.g., through increasing mitophagy) leads to lower mtDNA copy numbers and greater heteroplasmy variance. It is unclear whether a given treatment will have the sole effect of upregulating mitophagy: it seems likely that compensatory mechanisms (which we do not explicitly model, but may include retrograde signalling [Chae et al., 2013]) will engage to stabilise mtDNA copy number. However, such mechanisms would most straightforwardly be expected to act through increasing mtDNA proliferation, thus having the net effect of increasing mtDNA turnover. We have shown that such an increase in turnover also increases the heteroplasmy variance in a population. We therefore propose that upregulating mitophagy may be a fruitful pathway of investigation for increasing bottlenecking power, either as a standalone effect or due to the action of compensatory mechanisms it may invoke.

Speculatively, potential strategies to upregulate mitophagy may include the limited use of uncouplers to accelerate the mitophagy normally involved in quality control (de Vries et al., 2012); targetted chemical treatments with agents that have been identified as regulating mitophagy, including glutathione in yeast (Deffieu et al., 2009) and C₁₈-pyridium ceramide in human cancer cells (Sentelle et al., 2012); modulation of mitochondrial ultrastructure and dynamics to upregulate fission, intrinsically linked to the process of mitophagy (Youle and Narendra, 2010; Twig and Shirihai, 2011); or the use of existing drugs which have been found to modulate mitophagy, such as Efavirenz (Apostolova et al., 2011).

We have defined heteroplasmy by

To find statistics for this quantity we consider the Taylor expansion of a function f(X₁, X₂) of two random variables X₁, X₂ about a point (μ₁, μ₂), where ${\mu }_i=\mathbb{E}\left(X_i\right)$. We assume that the moments of X_i are well-defined and both have zero probability mass at X_i = 0. The Taylor expansion is:

where f_i denotes the derivative of f with respect to X_i. We truncate the expansion at first order for later algebraic simplicity, noting that even with this level of precision, the agreement between the resulting analysis and numerical simulation is excellent. Then

We note that $\mathbb{E}\left(X_i-{\mu }_i\right)=0$, so

Similarly,

and noting that $\mathbb{E}\left({\left(X_i-{\mu }_i\right)}^2\right)=\mathbb{V}\left(X_i\right)$ we obtain

where $\mathbb{C}\left(X_1,X_2\right)$ is the covariance of X₁ and X₂. If we now use $f\left(X_1,X_2\right)=\frac{X_1}{X_2}$, we have f₁ = X₂⁻¹, f₂ = −X₁X₂⁻²; so

If X₁ = M₂ and X₂ = M₁ + M₂, and M₁ and M₂ are independent (due to the lack of coupling between the mtDNA species), $\mathbb{C}\left(X_1,X_2\right)=\mathbb{V}\left(M_2\right)$, and so