Stochastic modelling, Bayesian inference, and new in vivo measurements elucidate the debated mtDNA bottleneck mechanism
Abstract
Dangerous damage to mitochondrial DNA (mtDNA) can be ameliorated during mammalian development through a highly debated mechanism called the mtDNA bottleneck. Uncertainty surrounding this process limits our ability to address inherited mtDNA diseases. We produce a new, physically motivated, generalisable theoretical model for mtDNA populations during development, allowing the first statistical comparison of proposed bottleneck mechanisms. Using approximate Bayesian computation and mouse data, we find most statistical support for a combination of binomial partitioning of mtDNAs at cell divisions and random mtDNA turnover, meaning that the debated exact magnitude of mtDNA copy number depletion is flexible. New experimental measurements from a wildderived mtDNA pairing in mice confirm the theoretical predictions of this model. We analytically solve a mathematical description of this mechanism, computing probabilities of mtDNA disease onset, efficacy of clinical sampling strategies, and effects of potential dynamic interventions, thus developing a quantitative and experimentallysupported stochastic theory of the bottleneck.
https://doi.org/10.7554/eLife.07464.001eLife digest
Mitochondria are structures that provide vital sources of energy in our cells. DNA contained within mitochondria encodes important mitochondrial machinery, and most human cells contain hundreds or thousands of mitochondrial DNA molecules in addition to the DNA that is stored in the nucleus. Mitochondrial DNA is inherited from mothers via the egg, and the details of this inheritance are poorly understood. This question is important because inherited mistakes in mitochondrial DNA can have detrimental consequences on health, with links to fatal diseases and many other conditions.
An unfertilised egg cell contains many copies of mitochondrial DNA molecules; some may have mutations and some may not. After fertilisation, the egg divides, the number of cells in the developing embryo increases, and the number of mitochondrial DNA molecules per cell changes. If the original egg cell contained defective mitochondrial DNA, some of these new cells end up containing more defective copies than others, leading to celltocell differences in the developing embryo. This potentially allows cells with the greatest number of defective mitochondria to be eliminated. The increase in this celltocell variability is called ‘bottlenecking’, and its mechanism remains highly debated.
Johnston et al. have now used tools from maths, statistics and new experiments to address this debate, in the light of several studies that measured the mitochondrial DNA content in developing mice. This approach allowed a new theoretical model of mitochondrial DNA during the growth of an organism to be produced, which encompasses a wide range of existing theories and allows them to be compared. This model starts from the viewpoint that the hundreds or thousands of mitochondrial DNA molecules in a cell can be thought of as a population undergoing random ‘birth’ and ‘death’, and it allows the first statistical comparison of the many proposed bottleneck mechanisms.
Johnston et al. find support for two ways that cells segregate mitochondria as they multiply, and show that the decrease in the number of mitochondrial DNA molecules during bottlenecking is flexible. This reconciles a debate amongst previous studies. These findings are confirmed using new experimental data from mice, which are genetically distinct from existing studies, illustrating the generality of the model's findings. Furthermore, an analytic mathematical description that describes in detail how bottlenecking might work is produced.
Finally, Johnston et al. provide examples using this new theoretical model to suggest therapeutic strategies for diseases caused by mitochondrial DNA mutations. Future work will need to test these suggestions, and link mathematical understanding of mitochondria with healthcare data.
https://doi.org/10.7554/eLife.07464.002Introduction
Mitochondria are vital energyproducing organelles within eukaryotic cells, possessing genomes (mitochondrial DNA, mtDNA) that replicate, degrade and develop mutations (Rand, 2001; Wallace and Chalkia, 2013). MtDNA mutations have been implicated in numerous pathologies including fatal inherited diseases and ageing (Lightowlers et al., 1997; Wallace, 1999; Poulton et al., 2009; Wallace and Chalkia, 2013). Combatting the buildup of mtDNA mutations is of paramount importance in ensuring an organism's survival. Substantial recent medical, experimental, and media attention has focused on methods to remove (Bacman et al., 2013) or prevent the inheritance of (Bredenoord et al., 2008; Poulton et al., 2009; Craven et al., 2010; Poulton et al., 2010; Burgstaller et al., 2015) mutated mtDNA in humans.
One means by which organisms may ameliorate the mtDNA damage that builds up through a lifetime is through a developmental process known as bottlenecking. Immediately after fertilisation, a single oocyte (which may contain >10^{5} individual mtDNAs) may have a nonzero mtDNA mutant load or heteroplasmy (the proportion of mutant mtDNA in the cell). As the number of cells in the developing organism increases, the intercellular population then acquires an associated heteroplasmy variance, that is, the variance in mutant load across the population of cells (Figure 1A), allowing removal of cells with high heteroplasmy and retention of cells with low heteroplasmy. Intense and sustained debate exists as to the mechanism by which this increase of heteroplasmy variance occurs. Several experimental results in mice suggest that, during development, the copy number of mtDNA per cell in the germ cell line drops dramatically to ∼10^{2}, reducing the effective population size of mitochondrial genomes (Cree et al., 2008; Wai et al., 2008). One postulated bottlenecking mechanism is that this low population size accelerates genetic drift and so increases the celltocell heteroplasmy variance (Bergstrom and Pritchard, 1998; Aiken et al., 2008; Cree et al., 2008; Wonnapinij et al., 2010), which was first observed to generally increase from primordial germ cells through primary oocytes to mature oocytes (Jenuth et al., 1996). However, independent experimental evidence (Wai et al., 2008) suggests that heteroplasmy variance increases negligibly during this copy number reduction, though this interpretation has been debated (Samuels et al., 2010). Wai et al. (2008) shows heteroplasmy variance rising during folliculogenesis, after the mtDNA copy number minimum has been passed. In yet another picture, supported by conflicting experimental results (Cao et al., 2007, 2009), heteroplasmy variance increases with a less pronounced decrease in mtDNA copy number (a minimum copy number >10^{3} in mice), solely through random effects associated with partitioning at cell divisions. Clearly a consensus on this important mechanism is yet to be reached.
Important existing theoretical work on modelling the bottleneck has assumed a particular underlying mechanism (Bergstrom and Pritchard, 1998; Wolff et al., 2011) or derived statistics of mtDNA populations (Chinnery and Samuels, 1999; Elson et al., 2001; Wonnapinij et al., 2008, 2010) without explicitly considering changing mtDNA population size, or the discrete nature of the mtDNA population: effects which may powerfully affect mtDNA statistics. To capture these effects it is necessary to employ a ‘bottomup’ physical description of mtDNA as populations of individual, discrete elements subject to replication and degradation, as in, for example, (Chinnery and Samuels, 1999) and (Capps et al., 2003). Exploring the bottleneck also requires explicitly modelling partitioning dynamics throughout a series of cell divisions, over which population size can change dramatically. While previous simulation work (Cree et al., 2008; Poovathingal et al., 2009) has taken such a philosophy with specific model assumptions, we are not aware of such a study allowing for the wide variety of replication and partitioning dynamics proposed in the literature; we further note that replicationdegradationpartitioning mtDNA models are yet to be fully described analytically. Nor is there a general quantitative framework under which different proposed bottleneck mechanisms can be statistically compared given extant data (although statistical analyses focusing on particular mechanisms and individual sets of experimental results have been used throughout the literature, for example, using a Bayesian approach under a particular bottleneck model to infer model bottleneck size [Marchington et al., 1998]). Combined developments in theory and inference are therefore required to make progress on this important question.
We remedy this situation by constructing a general model (features and parameters described in Figure 1) for the population dynamics of the bottleneck, able to describe the range of proposed mechanisms existing in the literature. Using experimental data on mtDNA statistics through development (Jenuth et al., 1996; Cao et al., 2007; Cree et al., 2008; Wai et al., 2008), we use approximate Bayesian computation (Beaumont et al., 2002; Toni et al., 2009; Sunnåker et al., 2013; Johnston, 2014) to rigorously explore the statistical support for each mechanism, showing that random mtDNA turnover coupled with binomial partitioning of mtDNAs at cell divisions is highly likely, and that the debated magnitude of mtDNA copy number reduction is somewhat flexible. Subsequently, we confirm the predictions of this model by performing new experimental measurements of heteroplasmy statistics in mice with an mtDNA admixture, including a wildderived haplotype, that is genetically distinct from previous studies. We then analytically solve the equations describing mtDNA population dynamics under this mechanism and show that these results allow us to investigate potential interventions to modulate the bottleneck (suggesting that upregulation of mtDNA degradation may increase the power of the bottleneck to avoid inherited disease; we discuss potential strategies for such an intervention) and yield quantitative results for clinical questions including the timescales and probabilities of disease onset, and the efficacy of strategies to sample heteroplasmy in clinical planning.
Results
A general mathematical model encompassing proposed bottlenecking mechanisms
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 birthdeathpartition (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 ‘bottomup’ 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 weaklyperfoming 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 nonselective 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 λ_{1}, ν_{1},…,λ_{6}, ν_{6}), 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, coarsegrained 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 suborganellar 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 highresolution 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 suborganellar mtDNA structure. We discuss the impact of mixed or fixed nucleoid content in Appendix 1.
A BDP model of mtDNA dynamics has most statistical support given experimental measurements
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 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 goodnessoffit 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 experimentallydetermined 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’).
Parameterisation and interpretation of the BDP model
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 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^{3}; 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^{−1} 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.
Experimental verification of the BDP model
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, wildderived 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 wildderived 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 nonzero, 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—source data 1
 https://doi.org/10.7554/eLife.07464.007
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 nonpathological 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 bestfit 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 AndersonDarling 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).
The BDP model is analytically tractable
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.
Mitochondrial turnover, degradation, and selective pressures exert quantifiable influence on heteroplasmy variance
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).
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, populationmaintaining 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 celllevel 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 λ_{1} ≠ λ_{2} (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’).
Probabilities of exceeding threshold heteroplasmy values
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_{1} wildtype and m_{2} 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.
Developmental sampling of embryonic 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 tissuespecific selective factors acting on the mtDNA population, this formalism provides a general means of rigorously inferring embryonic heteroplasmy through genetic diagnosis sampling.
Discussion
We have used a general stochastic model and approximate Bayesian computation with the available experimental data on developmental mtDNA dynamics to show that the bottleneck is most likely manifest through stochastic mtDNA dynamics and partitioning, with increased random turnover later during development, a mechanism which we can describe exactly and analytically (Figure 6). We emphasise that the bottomup construction of our model from physical first principles both increases the flexibility and generality of our model, allowing different mechanisms to be compared together, and providing information on mtDNA dynamics throughout development rather than estimating an overall effect. We note that even though our model cannot represent the full microscopic truth underlying the mtDNA bottleneck, its ability to recapitulate the wide range of extant experimental measurements suggest that its study may yield useful insights into the effects of different treatments and perturbations on the bottleneck.
A key debate in the literature has focussed on the magnitude of the bottleneck. Some studies (Aiken et al., 2008; Cree et al., 2008) have observed a depletion of mtDNA copy number during the bottleneck to minima around several hundred; other studies (Cao et al., 2007, 2009) have observed that mtDNA copy number remains >10^{3}. Our study shows that observed increases in heteroplasmy variance (Jenuth et al., 1996; Wai et al., 2008) can be achieved across this range of potential minimal mtDNA copy numbers, meaning that the muchdebated magnitude of mtDNA copy number reduction is not the sole critical feature of the bottleneck, in agreement with arguments from Cao et al. (2007, 2009); Wai et al. (2008). We find that the role of stochastic mtDNA dynamics can play a key role in determining heteroplasmy variance without additional mechanistic details, in keeping with approaches proposed by Cree et al. (2008). The mechanism with the most statistical support is thus consistent with aspects from all existing proposals in the literature.
We have shown that, of the models proposed in the literature, a BDP model, proposed after Cree et al. (2008) and compatible with an interpretation of Wai et al. (2008), is the individually most likely mechanism, and capable of producing experimentally observed heteroplasmy behaviour. We cannot, given current experimental evidence, discount hybrid mechanisms, where BDP dominates the population dynamics but small contributions from other mechanisms provide perturbations to this behaviour, and propose experiments to conclusively distinguish between these cases (see Appendix 1). As the expected statistics of mtDNA populations undergoing inheritance of heteroplasmic mtDNA clusters is very similar to those undergoing binomial partitioning of mtDNAs (see Appendix 1), the inheritance of heteroplasmic nucleoids (as opposed to individual mtDNAs) is not excluded by our findings, though other recent experimental evidence suggests that this situation may be unlikely (Poe et al., 2010; Jakobs and Wurm, 2014). We contend that the most likely situation may involve the partitioning of individual organelles, containing a mixture of homoplasmic nucleoids of characteristic size <2. Notably, this case (inheritance of heteroplasmic groups, likely with fluid structure due to mixing of organellar content and mitochondrial dynamics), gives rise to statistics which our binomial model reproduces (see Appendix 1).
As mentioned in the model description, it is likely that mitochondrial dynamics (fission and fusion of mitochondria) (Detmer and Chan, 2007) play a role in determining natural mtDNA turnover, and particularly mtDNA turnover in the presence of pathological mutations (Nunnari and Suomalainen, 2012), through the mechanism of mitochondrial quality control (Twig et al., 2008; Twig and Shirihai, 2011). Mitochondrial dynamics may also influence the elements of partitioning, through changes in the connectivity of the mitochondrial network. In our current model, these influences are coarsegrained into descriptions of the dynamic rates of mtDNA replication and degradation, and the characteristic elements that are partitioned at divisions. These physical parameters, as opposed to the more microscopic details of mitochondrial dynamics, are expected to be the key determinants of heteroplasmy statistics through development. Accounting for how these parameters, which summarize the relevant outputs of mitochondrial dynamics, connect to details of microscopic models of mitochondrial dynamics is an important future research direction to be addressed when more quantitative data is available.
The experimental data used to parameterise the first part of our study was taken from four studies in mice. Observation of similar dynamics in salmon (Wolff et al., 2011) points towards the bottleneck being a conserved mechanism in vertebrates. We also note that our results in mice are broadly consistent with findings from recent experiments in other organisms, suggesting that in primates and humans, heteroplasmy variance may increase at early developmental stages (Monnot et al., 2011; Lee et al., 2012), and that partitioning of mitochondria is binomial in HeLa cells (Johnston et al., 2012). As more studies become available on human mtDNA behaviour during development we will test our model's applicability and its clinical predictions. We note that the results of a recent study of human preimplantation sampling (Treff et al., 2012) found that earlier measurements provided strong predictive power of mean heteroplasmy, about which substantial variation was recorded in the offspring—both of which results are consistent with the application of our model to theoretical sampling considerations. In addition, recent observations that the m.3243A > G mutation in humans both increases mtDNA copy number during development (Monnot et al., 2013), and displays a less pronounced increase of heteroplasmy variance (Monnot et al., 2011) than other mutations, are consistent with the link between heteroplasmy variance and mtDNA copy number in our theory.
The combination of modern stochastic and statistical treatments that we have employed provides a generalisable and powerful way to recapitulate experimental data and rigorously deduce underlying biological mechanisms. We have used this combination to explore pertinent questions regarding the mtDNA bottleneck (and others have used a similar philosophy to numerically explore mtDNA point mutations [Poovathingal et al., 2009]): we hope to convince the reader that such methodology may be appropriate to explore other problems involving stochastic biological systems. We have used new experimental measurements to confirm our theoretical findings, illustrating the beneficial and powerful coupling of mathematical and experimental approaches to address competing hypotheses in the literature. Our detailed elucidation of the bottleneck allows us to propose further experimental methodology to address the current unknowns in our theory, including the specifics of mtDNA partitioning at cell division and the roles of selective differences between mtDNA types; importantly, we also propose a strategy to investigate our claim that our most supported model is compatible with the subsetreplication picture of mtDNA dynamics. We list these experiments in full in Appendix 1. Finally, we believe that the theoretical foundation for mtDNA dynamics that we have produced allows increased quantitative rigour in the predictions and strategies involved in mtDNA disease therapies, illustrated by the above application of our theory to problems in mtDNA sampling strategies, disease onset timescales, and interventions to increase the power of the bottleneck.
Materials and methods
General model for mtDNA dynamics
Request a detailed protocolOur ‘bottomup’ 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 nonwildtype 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 ‘wildtype’ mtDNAs in a cell as m_{1} and the number of ‘mutant’ mtDNAs as m_{2}. The heteroplasmy of a cell is then $h=\frac{{m}_{2}}{{m}_{1}+{m}_{2}}$, that is, the population proportion of mutant mtDNA.
MtDNA dynamics within a cell cycle
Request a detailed protocolIndividual 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.
MtDNA dynamics at cell divisions
Request a detailed protocolA parameter c (cluster size; a nonnegative 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.
Different dynamic phases through development
Request a detailed protocolThe 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).
Initial conditions
Request a detailed protocolThe initial conditions of our model involve an initial mtDNA copy number m_{0} (the total number of mtDNAs in the fertilised oocyte) and an initial heteroplasmy h_{0} (the fraction of these mtDNAs that are mutated).
Data acquisition
Request a detailed protocolWe 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]).
Simulation, model selection, and parametric inference
Request a detailed protocolWe 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^{3} 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^{3} simulated trajectories (all datapoints but one have n ≪ 10^{3}), 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 sumofsquares 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^{5} at all points throughout development, rejecting parameterisation that disobey this criterion. Metropolis MCMC was used to identify the bestfit 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 randomlyselected 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.
Creation of heteroplasmic mice
Request a detailed protocolHeteroplasmic 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 wildderived house mouse. Both mtDNA variants belong to the same subspecies, Mus musculus domesticus. For details on sequence divergence (see Burgstaller et al., 2014).
Isolation and lysis of oocytes
Request a detailed protocolMice were sacrificed at the indicated ages by cervical dislocation. Ovaries were extracted and immediately placed in cryobuffer 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 cryobuffer 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 96well PCR plates (Life Technologies, Austria) containing 10 μl of oocytelysis buffer (Lee et al., 2012) (50 mM TrisHCl, [p.H 8.5], 1 mM EDTA, 0.5% tween20 [Sigma–Aldrich, Austria] and 200 μg/ml Proteinase K [Macherey–Nagel, Germany]). Samples covered stages from primary oocytes of 3 dayold mice up to mature oocytes of 40 dayold 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 TrisEDTA buffer, pH 8.0 (Sigma–Aldrich, Austria). 3 μl were used per qPCR reaction.
Heteroplasmy quantification by Amplification Refractory Mutation System (ARMS)qPCR
Request a detailed protocolHeteroplasmy quantification was performed by ARMSqPCR, 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 realtime PCR experiments) guidelines (Bustin et al., 2009; Burgstaller et al., 2014). The proportion between HB derived and C57BL/6N mtDNA was determined by ARMSqPCR assays based on a SNP in mtrnr2 (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 wildderived <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.
Analytic model
Request a detailed protocolIn the BDP model, processes within a cell cycle constitute a birthdeath 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{\Sigma}}_{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 wildtype 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:
Selection
Request a detailed protocolThe 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_{0} is initial heteroplasmy and Δλ is the increase (or decrease, if negative) in replication rate of mutant over wildtype 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_{0} must lie between 0 and 1, and by the intuitive fact that heteroplasmy changes slow down as these limits are approached).
Threshold crossing
Request a detailed protocolThe 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 σ^{2} 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
Inference from heteroplasmy measurements
Request a detailed protocolGiven a sampled measurement heteroplasmy h_{m}, the probability $\mathbb{P}\left({h}_{0}{h}_{m}\right)$ that embryonic heteroplasmy is h_{0} 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 μ, σ^{2}, ζ_{1}, ζ_{2} are functions of h_{0}: μ, σ^{2} may be found numerically and the ζ values are analytically calculable (see Appendix 1).
Appendix 1
Data from experimental studies
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.
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).
Heteroplasmic and homoplasmic clusters
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 semipermanent 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 birthdeath 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_{0} mtDNAs, split into (1 − h)m_{0} wildtype and hm_{0} mutant molecules. These mtDNAs are clustered into m_{0}/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 celltocell heteroplasmy variance with time.
Appendix figure 1 shows the resultant heteroplasmy variance trajectories for different cases (with h_{0} = 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 longterm 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.
Parametric inference for bottlenecking dynamics
Our model is a function of the parameter set θ = {ν_{i}, λ_{i}, n_{i}, τ_{i}, S, α, T, c, h_{0}, m_{0}, δλ}. 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_{0} and initial copy number m_{0}; δλ 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^{−1}; S ∈ {0,1}; α ∈ [0.005,1]; T ∈ [0,100] day; c ∈ [0,100] day; h_{0} ∈ [0,1]; m_{0} ∈ [0,10^{6}]. The following values are fixed from experimental studies (Lawson and Hage, 1994; Drost and Lee, 1995): n_{1} = 29; n_{2} = 7; τ_{1} = 7 hr; τ_{2} = 16 hr. λ_{6} = 0 hr^{−1} is fixed to avoid mtDNA proliferation after development; h_{0} = 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{\u03f5}\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 sumofsquares 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_{1} = 10^{3}. 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 sumofsquares 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_{1} over orders of magnitude from 10^{2}–10^{4} (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_{1} = 3 × 10^{3} 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{\u03f5}$. We ran 10^{6} 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^{6} 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.
Posteriors for all variables and datasets
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.
Experimental measurements
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:
Co2f: GCCAATAGAACTTCCAATCCGTATAT,
Co2r: TGGTCGGTTTGATGTTACTGTTG,
Co2FAM: CTGATGCCATCCCAGGCCGACTAABHQ1 (Amplicon length: 136 bp);
Co3f: TCTTATATGGCCTACCCATTCCAA,
Co3r: GGAAAACAATTATTAGTGTGTGATCATG,
Co3FAM: TTGGTCTACAAGACGCCACATCCCCTBHQ1 (Amplicon length: 103 bp).
ARMSassays:
16SrRNA2340(3)Gf: AATCAACATATCTTATTGACCaAG (haplotype C57BL/6N),
16SrRNA2340(3)Af: AATCAACATATCTTATTGACCgAA (haplotype HB);
16SrRNA2458r: CAC CAT TGG GAT GTC CTG ATC,
16SrRNAFAM: FAMCAA TTA GGG TTT ACG ACC TCG ATG TTBHQ1. (Amplicon length: 142 bp).
Every qPCR run consisted of one consensus and an ARMS assay.
Mastermixes for triplicate qPCR reactions contained 1× buffer B (Solis BioDyne, Estonia); 4.5 MgCl_{2} for the ARMS and the Co3 consensus assays, and 3.5 mM MgCl_{2} 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 mastermix and 3 μl DNA were transferred in triplicates to 384well PCR plates (Life Technologies, Austria) using the automated pipetting system epMotion 5075TMX (Eppendorf, Germany). Amplification was performed on the ViiA 7 RealTime 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 wildderived 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 Yintercept = 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 nonanalysed 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 TrisEDTA buffer (Sigma–Aldrich, Austria), pH 8.0.
Bottlenecking mechanisms and further experimental elucidation
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.
Other proposed mechanisms
Random partitioning of homoplasmic mtDNA clusters
Cao et al. (2007) suggests a less powerful depletion of mtDNA copy number during early development than assumed by other studies, with heteroplasmy variance increase instead being explained by the partitioning of clusters of mtDNA at cell divisions. However, the time period over which Wai et al. (2008) and Jenuth et al. (1996) observe increasing heteroplasmy variance corresponds to a situation in which germ line cells are largely quiescent, immediately suggesting that partitioning at cell divisions cannot explain increasing variance (as cell divisions do not occur). Furthermore, results from our model suggest that, unless these clusters are very small, this mechanism would immediately lead to a rather higher and sharper increase in heteroplasmy variance than observed.
Replication of a specific subset of mtDNAs during folliculogenesis
Wai et al. (2008) proposes a mechanism in which only a subset of mtDNAs replicate during folliculogenesis. There are several specific dynamic schemes by which this mechanism could be manifest. The first that we consider involves the following scenario: at some point during development, around the start of folliculogenesis, a specific subset of mtDNAs in each cell is ‘marked’ as able to replicate (we next consider the case in which this subset is more plastic with time). In this case, the effect of ‘switching off’ replication of a subset of mtDNAs depends on the balance of replication and degradation rates of the mtDNA population:
Low replication, low degradation. In this case, the population stays largely static; the switching off of replication has little effect, and the heteroplasmy variance cannot increase to the levels observed in experiment.
Low replication, high degradation. In this case, the high degradation rate ensures that the nonreplicating mtDNAs are removed from the cell, providing a ‘bottleneck’ as only the replicating mtDNAs remain. However, this regime yields a transient period of mtDNA copy number depletion, while the nonreplicating mtDNAs are degrading but the (small) population of replicating agents remains low. This copy number depletion is not observed.
High replication, high degradation. In this case, nonreplicating mtDNAs are removed and replicating mtDNAs are capable of fast enough replication to survive the transient drop in copy number. However, the rates associated with this mechanism are necessarily high enough such that the increase in heteroplasmy variance is very sharp, notably more so than the smooth increase with time observed in experiment (see, e.g., Figure 2 in the Main text).
Combined subset replication, and/or heteroplasmic cluster inheritance, with and random dynamics. Our approach does not provide support against a model combining our inferred mechanism (increased random turnover of mtDNAs) with some other dynamic schemes, namely (a) that in which only a subset of mtDNAs may replicate during folliculogenesis, and/or (b) where heteroplasmic mtDNA clusters, rather than individual mtDNAs, are the units of inheritance. A combination with (a) would allow the reduction of the key parameters associated with each: so the rate of random turnover could be lower, and the proportion of replicating genomes larger, than in the case of the pure incarnations of those respective mechanisms. This scheme may thus provide a viable alternative—however, it requires an introduction of two coupled mechanisms, which experimental data currently cannot disambiguate. For this reason and for parsimony, we report the case where random dynamics alone are responsible, and below suggest experimental protocols to further elucidate this possible link or its absence. A combination with (b) is possible and cannot be discounted using the available data, as the trajectories of heteroplasmy variance under (b) and under binomial inheritance are the same. We propose observations of mitochondrial ultrastructure and mtDNA localisation during development to resolve this remaining mechanistic question.
Observation of a subset of replicating genomes
Wai et al. (2008) performs BrU labelling to observe the proportion of mitochondria replicating in primary oocytes between P14 (21–25 dpc on our time axis). The observations contained therein (Appendix figure 2 in Wai et al., 2008) show a small subset of BrUlabelled 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^{−1} 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^{−1} 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.
Experimental elucidation
In Table 3 we list several classes of potential experimental protocols that would assist in further elucidation of the bottlenecking mechanism and our predictions. Potentially useful results include further characterisation of the microscopic detail underlying mtDNA dynamics during development, confirmation of our random turnover model, assessing degree to which heteroplasmy modulates copy number dynamics and exploring our predictions relating mitophagy and bottlenecking power.
Mitophagy regulation
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_{18}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).
Heteroplasmy statistics
We have defined heteroplasmy by
To find statistics for this quantity we consider the Taylor expansion of a function f(X_{1}, X_{2}) of two random variables X_{1}, X_{2} about a point (μ_{1}, μ_{2}), where ${\mu}_{i}=\mathbb{E}\left({X}_{i}\right)$. We assume that the moments of X_{i} are welldefined 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_{1} and X_{2}. If we now use $f\left({X}_{1},{X}_{2}\right)=\frac{{X}_{1}}{{X}_{2}}$, we have f_{1} = X_{2}^{−1}, f_{2} = −X_{1}X_{2}^{−2}; so
If X_{1} = M_{2} and X_{2} = M_{1} + M_{2}, and M_{1} and M_{2} 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
Derivation of analytic results for binomial model
Generating function within a cell cycle
To make analytic progress describing the mitochondrial content of quiescent cells, and within a single cell cycle of dividing cells, we use a birth and death model to describe mitochondrial evolution. Without cell divisions, the dynamics of a population of replicating and degrading entities is given by the master equation
with P(m) the probability of observing the system with a copy number m at time t, and m_{0} the initial copy number. The corresponding generating function, using the transformation $G\left(z,t\right)={\displaystyle \sum}_{m}{z}^{m}P\left(m,t\right)$, obeys
which is straightforwardly solved by
where the 0 subscript signifies that no divisions have occurred, and we have specifically labelled the base of G_{0} as g_{0} for later convenience.
Generating function over cell divisions
We now consider a system undergoing cell divisions. Now, we have a population of organelles with time evolution described by a generating function $G={\left[g\right]}^{{m}_{0}}$ and subject to binomial partitioning at cell division. The probability distribution of m after a single cell division is:
where m_{i,a}, m_{i,b} mean respectively the number of individuals after and before the ith cell division, and the subscript in P_{0} denotes the fact that this function refers to time evolution within a cell cycle (with no division). The sum takes into account all possible configurations of the system up to the cell division then all possible configurations afterwards, with weighting according to a binomial partitioning. This line of reasoning can straightforwardly be extended to n cell divisions (Rausenberger and Kollmann, 2008):
where Φ_{i} is a ‘probability propagator’ of the form
and m_{n+1,a} ≡ m_{0}. For clarity, we introduce the nomenclature:
Now consider the generating function of P_{n}:
where we have used the identity ${{{\displaystyle \sum}}^{\text{}}}_{a=0}^{b}{x}^{a}\left(\begin{array}{c}b\\ a\end{array}\right){2}^{b}\equiv {\left(\frac{1}{2}+\frac{x}{2}\right)}^{b}$ and changed variables $z\prime =\frac{1}{2}+\frac{{g}_{0}\left(z,t\right)}{2}$. Comparing Equations 29, 32 and following this process by induction we can see that the overall generating function is ${G}_{n}={h}_{0}^{{m}_{0}}$, where h is the solution to the recursive system
h_{i} is of the form $\frac{a{h}_{i+1}+b}{c{h}_{i+1}+d}$ (from Equation 21). This expression takes the form of a Riccati difference equation and can be solved exactly after Brand (1955). The solution is straightforward but algebraically lengthy, and we defer presentation of the full procedure to a future technical publication (Johnston and Jones, 2015). The overall solution is:
where the C subscript denotes cycling cells, and
Generating function for different phases
We now consider how to extend this reasoning to the overall bottlenecking process, which in general may involve several phases of quiescent and cycling dynamics with different kinetic parameters. We begin with the generating function bases g_{i}(z, t) for each regime i. For consistency with the above approach, we label phases starting from a zero index, so the first phase corresponds to i = 0, and we use i_{max} to denote the label of the final phase. Then we use
using induction over the different phases in the way we used induction over different cell cycles above. Here we consider the changeover between regimes by using the generating function at the start of the incoming phase.
The appropriate generating function bases for quiescent (Equation 21) and cycling (Equation 35) cells can be written as
with coefficients
using, as before, l = e^{(λ − ν)τ} and l′ = e^{(λ − ν)t}. Note that the cycling coefficients reduce to the quiescent coefficients when n → 0 and τ → 0. The values of the appropriate A, B, C, D coefficients for a given dynamic phase thus follow straightforwardly from the kinetic parameters of that phase, with the appropriate choice between quiescent and cycling parameters being made.
If we now label these coefficients with a subscript denoting the appropriate phase of bottlenecking, so that, for example, A_{i} is Equation 47 with λ_{i}, ν_{i}, n_{i} replacing λ, ν, n, we can write:
Following this recursion for n phases of bottlenecking and simplifying the resultant multilayer fraction gives rise to the solution
where
from which ${G}_{overall}={g}_{overall}^{{m}_{0}}$ follows straightforwardly. The following results will be of assistance:
As Equations 43–46 can be thought of as special cases of Equations 47–50, we combine Equations 47–50 into Equation 55, and, simplifying, we find the following relations:
we then immediately obtain
leaving us only with the problem of calculating the expression $\left(B\prime C\prime \left({m}_{0}+1\right)+A\prime \left(2C\prime +D\prime \left(1{m}_{0}\right)\right)\right)$ in the variance calculation. We were not able to dramatically simplify this expression and so, for clarity, write:
which gives us:
We note that Φ is just a notational simplification and is straightforwardly calculable by inserting Equations 47–50 into Equation 55 then computing Equation 64.
Constant population size
For generality, we consider enforcing a constant population size in postmitotic cells (not undergoing divisions). This process involves setting λ = ν, so the net gain in mtDNA is zero. If we write λ = ν + ϵ and take the limit ϵ → 0, Equation 21 becomes
To enforce a constant mean population size in mitotic cells, it is necessary to balance the expected loss of mtDNA through repeated divisions with an expected increase during the cell cycle. This balance can be accomplished by setting $\lambda =\nu +\frac{\text{ln}2}{\tau}$. Writing $\lambda =\nu +\frac{\text{ln}2}{\tau}+\mathit{\u03f5}$ and taking the ϵ → 0 limit we obtain
In both these cases, the same approach as above can be used to derive moments of the resulting probability distributions.
Explicit distributions
The probability of observing exactly m mtDNAs of a given type can be found from the generating function with
We can use Leibniz' rule on a generating function of form $G={\left(\frac{A\prime z+B\prime}{C\prime z+D\prime}\right)}^{{m}_{0}}$ by setting $f\equiv {\left(A\prime z+B\prime \right)}^{{m}_{0}}$, $g\equiv {\left(C\prime z+D\prime \right)}^{{m}_{0}}$ and writing
Enforcing z = 0 and rewriting in terms of a hypergeometric function gives
The distribution of heteroplasmy is then given by
where $\mathcal{I}\left(h\prime ,h\right)$ is an indicator function returning 1 if h′ = h and 0 otherwise. Computing the probability of observing a given heteroplasmy thus involves a sum, over all mtDNA states that correspond to that heteroplasmy, of the probability of that state.
The evaluation of hypergeometric functions is more computationally demanding than that of more common mathematical functions, and the infinite sums at first glance seem intractable. However, in practise and using parameterisations from our inferential approach, vanishingly little probability density exists at m_{1}, m_{2} > 5 × 10^{5}, corresponding to the biological observation that mtDNA copy number is very unlikely to exceed this value. Dynamic programming then allows these sums to be performed straightforwardly.
Finally, the computation of $\mathbb{P}\left(m=0,t\right)$ is important in our analysis of the characterisation of key distributions using the first two moments (see below), where it appears as $\mathbb{P}\left({m}_{2}=0,t\right)$, the probability of wildtype fixation. This is relatively straightforward to address analytically as when m = 0, Equation 68 reduces to $\mathbb{P}\left(0,t\right)={G}_{overall}{}_{z=0}$, which in the notation above is simply:
where we introduce the notation ζ for fixation probability for later brevity. We could not dramatically simplify the full expression so we leave it in this form and note that it can be readily calculated (as above) by inserting Equations 47–50 into Equation 55 then computing Equation 73.
Multiple species and heteroplasmy
The heteroplasmy h = m_{2}/(m_{1} + m_{2}) is straightforwardly addressable by considering the above solutions for m_{1} and m_{2}. We can also consider a more general case, in which we have four species of mtDNA in our model: wildtype reproducing (m_{1}), mutant reproducing (m_{2}), wildtype sterile (m_{3}) and mutant sterile (m_{4}). We assume that these species evolve in an uncoupled way with time. The parameter h_{0}, initial heteroplasmy, determines the initial proportion of mutant genomes: ${h}_{0}=\frac{{m}_{20}+{m}_{40}}{{m}_{0}}$, where m_{0} = m_{10} + m_{20} + m_{30} + m_{40} is the total initial copy number of mtDNA. The parameter α determines the proportion of genomes capable of reproducing: $\alpha =\frac{{m}_{10}}{{m}_{10}+{m}_{30}}=\frac{{m}_{20}}{{m}_{20}+{m}_{40}}$. We compute the time trajectories for all m_{i} then calculate heteroplasmy by setting M_{1} = m_{1} + m_{3}, M_{2} = m_{2} + m_{4}, respectively the total numbers of wildtype and mutant mtDNAs, and using Equations 15, 16, where all means and variances are straightforwardly extracted from the above analysis.
Characterisation of distributions of important quantities with moments
We are interested in the probability with which heteroplasmy h exceeds a certain threshold value h^{*}. This probability can be computed using Equation 72 above, but the large sums of hypergeometric functions suggest that a simpler approximation of the heteroplasmy distribution may be desirable, both for computational simplicity and intuitive interpretability. We here explore how well distributions of copy number and, importantly, heteroplasmy are characterised by quantities that are easily obtained from our analytic approaches without large summations: specifically, loworder moments $\mathbb{E}\left(m\right),\mathbb{V}\left(m\right)$, and fixation probabilities $\mathbb{P}\left(m=0\right)$.
For moderate initial heteroplasmy 0.7 > h_{0} > 0.3, all distributions are well matched by the Normal distributions computed using the first two moments $\mathbb{E}\left(m\right)$ and $\mathbb{V}\left(m\right)$. This match begins to fail as initial heteroplasmy decreases or increases to the extent where fixation of one mtDNA type becomes likely. The resultant nonnegligible probability density at h = 0 and/or h = 1 represents a truncation point which forces skew on the distributions (particularly $\mathbb{P}\left(h\right)$) and weakens the Normal approximation.
We can make progress by considering $\mathbb{P}\left(h\right)$ to be a weighted sum of a truncated Normal distribution $\mathcal{N}\prime \left(\mu ,{\sigma}^{2}\right)$ (truncated at 0, 1; and with currently unknown parameters μ, σ) and two δfunctions at h = 0 and h = 1 representing the fixation probability of wildtype and mutant mtDNA respectively. If we write ${\mathbb{P}}_{\mathcal{N}\prime}\left(h\right)$ for the probability density at h of such a truncated Normal distribution, we have:
where ${\zeta}_{1}=\mathbb{P}\left({m}_{2}=0,t\right)$ is the fixation probability of the wildtype and ${\zeta}_{2}=\mathbb{P}\left({m}_{1}=0,t\right)$ is the fixation probability of the mutant, expressions for which were computed previously in Equation 73. Knowledge of the parameters μ, σ that describe the truncated Normal part of this distribution will then provide us with a better estimate of $\mathbb{P}\left(h\right)$.
We can use the relations $\mathbb{E}\left(h\right)={{\displaystyle \int}}^{\text{}}h\mathbb{P}\left(h\right)dh$ and $\mathbb{V}\left(h\right)={{\displaystyle \int}}^{\text{}}{h}^{2}\mathbb{P}\left(h\right)dh\mathbb{E}{\left(h\right)}^{2}$. As δ(h) provides a nonzero contribution to these integrals only when h = 0, the contribution from this part of $\mathbb{P}\left(h\right)$ is always zero; then,
where $\mathbb{E}\left(\mathcal{N}\prime \right)$, $\mathbb{V}\left(\mathcal{N}\prime \right)$ are respectively the mean and variance of the truncated Normal distribution, and in the final line we have used the fact that $\mathbb{V}\left(\mathcal{N}\prime \right)={{\displaystyle \int}}^{\text{}}{h}^{2}{\mathbb{P}}_{\mathcal{N}\prime}\left(h\right)dh\mathbb{E}{\left(\mathcal{N}\prime \right)}^{2}$.
Results are known (Greene, 2003) for moments of the truncated Normal distribution:
where, in our case (with truncations at h = 0 and h = 1) α_{1} = −μ/σ, α_{2} = (1 − μ)/σ and $f\left(x\right)={\left(\sqrt{2\pi}\right)}^{1}\text{exp}\left({x}^{2}/2\right)$ and $F\left(x\right)=\frac{1}{2}\left(1+\text{erf}\left(x/\sqrt{2}\right)\right)$ are respectively the p.d.f. and c.d.f. of the standard Normal distribution. Given these expressions, we wish to invert these Equations 76, 78 to find μ and σ, the parameters underlying the truncated Normal distribution, given $\mathbb{E}\left(h\right)$, $\mathbb{V}\left(h\right)$ and ${\zeta}_{\mathrm{1,2}}=\mathbb{P}\left({m}_{\mathrm{2,1}}=0\right)$, which we can compute (see below). We have not been able to find an analytic solution for these equations; however, numerically solving these equations is computationally far cheaper than performing the numeric simulations required to better characterise the real distribution. We then obtain an expression for $\mathbb{P}\left(h\right)$, which well matches the exact distribution derived using Equation 72 (see Appendix figure 4).
Threshold crossing
The probability of crossing a threshold heteroplasmy h^{*} with time is simply given by the probability density in the region h > h^{*}. We can then use the result
for threshold crossing, which follows straightforwardly from considering the integrated density of the model distribution (Equation 74) of h above h^{*}, with parts from the error function representing the definite integral of the truncated Normal part of the distribution, with additional terms from wildtype fixation (if h^{*} ≠ 0) and mutant fixation (if h^{*} ≠ 1).
Inferring embryonic heteroplasmy
The probability that a sample measurement h_{m} came from an embryo with heteroplasmy h_{0} can be found from Bayes' Theorem:
We assume a uniform prior distribution $\mathbb{P}\left({h}_{0}\right)=\rho $ on embryonic heteroplasmy (though this can be straightforwardly generalised). $\mathbb{P}\left({h}_{m}\right)$ is given by the integral over all possible embryonic heteroplasmies of making observation h_{m}, so we obtain
where the μ, σ^{2} moments characterising the truncated Normal distribution are found numerically as above (for each ${h}_{0}\prime $ value in the integrand, which is performed numerically); and ζ_{1}, ζ_{2} are also functions of h_{0}.
References

Variations in mouse mitochondrial DNA copy number from fertilization to birth are associated with oxidative stressReproductive Biomedicine Online 17:806–813.https://doi.org/10.1016/S14726483(10)604099

Germline bottlenecks and the evolutionary maintenance of mitochondrial genomesGenetics 149:2135–2146.

Mitochondrial DNA nucleoid structureBiochimica et Biophysica Acta 1819:914–920.https://doi.org/10.1016/j.bbagrm.2011.11.005

A sequence defined by a difference equationThe American Mathematical Monthly 62:489–492.https://doi.org/10.2307/2307362

Ooplasmic and nuclear transfer to prevent mitochondrial DNA disorders: conceptual and normative issuesHuman Reproduction Update 14:669–678.https://doi.org/10.1093/humupd/dmn035

Mitochondrial DNA disease and developmental implications for reproductive strategiesMolecular Human Reproduction 21:11–22.https://doi.org/10.1093/molehr/gau090

A model of the nuclear control of mitochondrial DNA replicationJournal of Theoretical Biology 221:565–583.https://doi.org/10.1006/jtbi.2003.3207

Relaxed replication of mtDNA: a model with implications for the expression of diseaseThe American Journal of Human Genetics 64:1158–1165.https://doi.org/10.1086/302311

Glutathione participates in the regulation of mitophagy in yeastJournal of Biological Chemistry 284:14828–14837.https://doi.org/10.1074/jbc.M109.005181

Functions and dysfunctions of mitochondrial dynamicsNature Reviews Molecular Cell Biology 8:870–879.https://doi.org/10.1038/nrm2275

Biological basis of germline mutation: comparisons of spontaneous germline mutation rates among Drosophila, mouse, and humanEnvironmental and Molecular Mutagenesis 25:48–64.https://doi.org/10.1002/em.2850250609

Random intracellular drift explains the clonal expansion of mitochondrial DNA mutations with ageThe American Journal of Human Genetics 68:802–806.https://doi.org/10.1086/318801

Mitochondrial nucleoids maintain genetic autonomy but allow for functional complementationThe Journal of Cell Biology 181:1117–1128.https://doi.org/10.1083/jcb.200712101

Exact stochastic simulation of coupled chemical reactionsThe Journal of Physical Chemistry 81:2340–2361.https://doi.org/10.1021/j100540a008

Integration of cellular bioenergetics with mitochondrial quality control and autophagyBiological Chemistry 393:1485–1512.https://doi.org/10.1515/hsz20120198

Superresolution microscopy of mitochondriaCurrent Opinion in Chemical Biology 20:9–15.https://doi.org/10.1016/j.cbpa.2014.03.019

Efficient parametric inference for stochastic biological systems with measured variabilityStatistical Applications in Genetics and Molecular Biology 13:379–390.https://doi.org/10.1515/sagmb20130061

Closedform stochastic solutions for nonequilibrium dynamics and inheritance of cellular components over many cell divisionsProceedings of the Royal Society A In press.

Mitochondrial variability as a source of extrinsic cellular noisePLOS Computational Biology 8:e1002416.https://doi.org/10.1371/journal.pcbi.1002416

Solution of a process of random genetic drift with a continuous modelProceedings of the National Academy of Sciences of USA 41:144–150.https://doi.org/10.1073/pnas.41.3.144

mtDNA makes a uturn for the mitochondrial nucleoidTrends in Cell Biology 23:457–463.https://doi.org/10.1016/j.tcb.2013.04.009

Clonal analysis of the origin of primordial germ cells in the mouseGermline Development 165:68–84.

Mammalian mitochondrial genetics: heredity, heteroplasmy and diseaseTrends in Genetics 13:450–455.https://doi.org/10.1016/S01689525(97)012663

Evidence from human oocytes for a genetic bottleneck in an mtDNA diseaseThe American Journal of Human Genetics 63:769–775.https://doi.org/10.1086/302009

Markov chain Monte Carlo without likelihoodsProceedings of the National Academy of Sciences of USA 100:15324–15328.https://doi.org/10.1073/pnas.0306899100

Mutation dependance of the mitochondrial DNA copy number in the first stages of human embryogenesisHuman Molecular Genetics 22:1867–1872.https://doi.org/10.1093/hmg/ddt040

Detection of heteroplasmy in individual mitochondrial particlesAnalytical and Bioanalytical Chemistry 397:3397–3407.https://doi.org/10.1007/s0021601037513

Stochastic drift in mitochondrial DNA point mutations: a novel perspective ex silicoPLOS Computational Biology 5:e1000572.https://doi.org/10.1371/journal.pcbi.1000572

The units of selection of mitochondrial DNAAnnual Review of Ecology and Systematics 32:415–448.https://doi.org/10.1146/annurev.ecolsys.32.081501.114109

Quantifying origins of celltocell variations in gene expressionBiophysical Journal 95:4523–4528.https://doi.org/10.1529/biophysj.107.127035

Reassessing evidence for a postnatal mitochondrial genetic bottleneckNature Genetics 42:471–472.https://doi.org/10.1038/ng0610471

Ceramide targets autophagosomes to mitochondria and induces lethal mitophagyNature Chemical Biology 8:831–838.https://doi.org/10.1038/nchembio.1059

Approximate Bayesian computationPLOS Computational Biology 9:e1002803.https://doi.org/10.1371/journal.pcbi.1002803

Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systemsJournal of the Royal Society Interface 6:187–202.https://doi.org/10.1098/rsif.2008.0172

Blastocyst preimplantation genetic diagnosis (PGD) of a mitochondrial DNA disorderFertility and Sterility 98:1236–1240.https://doi.org/10.1016/j.fertnstert.2012.07.1119

Mitochondrial fusion, fission and autophagy as a quality control axis: the bioenergetic viewBiochimica et Biophysica Acta 1777:1092–1097.https://doi.org/10.1016/j.bbabio.2008.05.001

The interplay between mitochondrial dynamics and mitophagyAntioxidants & Redox Signaling 14:1939–1951.https://doi.org/10.1089/ars.2010.3779

Mitochondrial diseases in man and mouseScience 283:1482–1488.https://doi.org/10.1126/science.283.5407.1482

Mitochondrial DNA genetics and the heteroplasmy conundrum in evolution and diseaseCold Spring Harbor Perspectives in Biology 5:a021220.https://doi.org/10.1101/cshperspect.a021220

The distribution of mitochondrial DNA heteroplasmy due to random genetic driftThe American Journal of Human Genetics 83:582–593.https://doi.org/10.1016/j.ajhg.2008.10.007

Mechanisms of mitophagyNature Reviews Molecular Cell Biology 12:9–14.https://doi.org/10.1038/nrm3028

Mitochondrial fission, fusion, and stressScience 337:1062–1065.https://doi.org/10.1126/science.1219855
Decision letter

Jodi NunnariReviewing Editor; University of California, Davis, United States
eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.
Thank you for resubmitting your work entitled “Stochastic modelling, Bayesian inference, and new in vivo measurements elucidate the debated mtDNA bottleneck mechanism” for further consideration at eLife. Your revised article has been favorably evaluated by Aviv Regev (Senior editor) and a member of the Board of Reviewing Editors along with two reviewers. The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance as detailed in the reviewers' comments below.
In particular, as pointed out in the previous version, claims of novelty in several areas are overstated and require the appropriate revision for accuracy. The reviewers' concerns about data variance in the experimental data also need to be addressed (especially in the 25 day old data) as well as the concern that there is a very large effect with only small alterations in mtDNA degradation.
Reviewer #1:
This is a resubmission of a pure modeling paper, with the addition of new experimental data which is presented as supporting the main conclusions of the modeling. The authors have carried out a very thorough theoretical and analytical analysis of the proposed mtDNA genetic bottleneck, and conclude that the existing stance taken by Cree and Wai model (b) is a sensible interpretation of the experimental data. As such, the current paper is reassuringly confirmatory. However, I still have concerns, in part relating to the way the paper is written, and in part relating to the new data, which raises additional questions.
Introduction. A statement is made: ‘No study yet exists combining an analytic ‘bottomup’ physical description of the bottleneck with mtDNAs as individual, discrete elements subject to a possible variety of replication and partitioning dynamics throughout an explicitly modelled series of cell divisions’. This is not correct, and I made this point in my earlier review of the paper. Their reference (Cree et al., 2008) includes such a ‘bottom up’ model. Although the compartments did not vary in that model, assumptions were made about the partitioning which are in accordance with the conclusions of this current paper.
Likewise, the statement ‘theory describing the behaviour of mtDNA populations throughout development, as opposed to estimating an overall effect of the bottleneck, is also absent’ is also an overstatement designed to justify the current work. There are several such theories that the authors refer to in their Introduction.
For the new experimental data, they use a line previously reported in their Cell Reports publication. A key message of this publication was that different mtDNA haplotypes behave differently, that selection occurs, and that this is influenced by the nuclear genetic background. All of these effects will influence the heteroplasmy variance. Although the authors touch on selection, they do not explore this thoroughly in the context of their previous work. How will this impact on their interpretation of their bottleneck models?
For the new experimental data, the heteroplasmy measurements are reported to a remarkable degree of accuracy. Do the authors really think that they can reliably report down to 0.01% as indicated in the supplementary data?
Previous theoretical work has shown that ∼50 heteroplasmy measurements are required to reliably compute heteroplasmy variance. Several of the published datasets do not achieve this, and thus the reported variance measurements are likely to be inaccurate. The same may also apply to the new experimental dataset.
Figure 4C raises concerns. For the ∼25 day old data, 6 data points (from a total of 15) are way outside the confidence intervals predicted by the various models. This implies that the new data set is not supporting the current or previous models.
Finally, given the likely readership, it is important to keep reminding readers what is prediction, and what is actual measurement. The section on turnover illustrates this nicely. No turnover measurements have been made; indeed, they would be incredibly difficult to do.
Reviewer #2:
The authors have revised their previous submission to now include new experimental data on the mtDNA bottleneck, as requested in an earlier round of review.
There is clearly a great deal of variation in the new data for the normalized variance, but the new data does fit best with the authors' favoured birthdeathpartition model. The inclusion of this extra data, at the time points which are most discriminatory with regard to the various models, is therefore helpful, though not absolutely definitive. In particular, it looks as if the spread of the normalised variance itself has interesting dynamics, starting out low, increasing to a maximum at around 25 dpc before diminishing again. Is there any way in which the BDP model can capture these dynamics? Also it was not clear to me why the BDP curve in Figure 4A differs from the HB theory mean plot in Figure 4C.
Finally, I remain concerned that a 2% alteration in mtDNA degradation can have such a large effect, see Figure 5B. The authors must address this issue.
[Editors’ note: a previous version of this study was rejected after peer review, but the authors submitted for reconsideration. The previous decision letter after peer review is shown below.]
Thank you for choosing to send your work entitled “Stochastic modelling and Bayesian inference elucidate the debated mechanism of the mtDNA bottleneck” for consideration at eLife. Your full submission has been evaluated by Aviv Regev (Senior editor), a Reviewing editor and three peer reviewers, and the decision was reached after discussions between the reviewers. We regret to inform you that your work will not be considered further for publication.
Although there was enthusiasm for the quantitative approach to the mitochondrial DNA bottleneck, the reviewers felt that experiments directed at testing the predictions of the modeling are required for the manuscript to significantly advance the field. Therefore, at this stage, we are returning the manuscript to you with an encouragement to submit with experimental data. We hope the reviewers' comments are helpful.
Reviewer #1:
The manuscript from Johnston et al. is an extremely thorough and wellthought out attempt to inject rigorous, quantitative thinking into the mechanism of the mtDNA bottleneck. It would appear that prior attempts in this direction have been limited in scope and inconclusive. The strength of the authors' approach is that it is able to unify previous treatments of the phenomenon into a single theoretical framework. From this position, it is possible to assess the plausibility of the various mechanisms within a Bayesian statistical analysis. Using prior data, the authors conclude that a birthdeathpartition model is most probable. This model is analytically tractable which is also advantageous. These are powerful conclusions, and it would seem absolutely certain that progress in this field is going to require an analysis of this sort, due to the inherent complexity (and stochasticity) of the problem.
I do however have some concerns. Principally, the analysis is based entirely on previously published data, which is not all internally consistent. In particular, looking at Figure 2A for the variances, there doesn't seem to be very good agreement between the different datasets (especially comparing the Jenuth data with the rest). I appreciate that the authors' statistical analysis takes into account the experimental uncertainty, but the divergence between the different experimental data sets found in the different studies suggests that the effective experimental error bars may be much bigger than estimated. The ability to distinguish between the models via Bayesian model comparison may then be diminished. This point is important, as there is relatively little difference between the models with respect to the mean values (whose experimental values are also well clustered between the various data sets). Hence it does seem that the variances are the key to distinguishing between the models. Also the statistical inference is based around a residual sum of squares difference. However, the way this was implemented (using logs of the copy number measurements, for example) seemed a bit arbitrary. Were other possibilities investigated to ensure consistency?
More generally, I am concerned that the modelling is used primarily as a tool to fit existing data and thereby determine which model is most likely on the basis of previously published data. It would be more powerful if the authors could actually test the novel predictions of their modelling in new experiments. This would elevate the study to a higher level. Their model most definitely is predictive (Table 2), it's just that they ought to consider teaming up with an experimental group and prove these predictions right!
One other issue concerns the sensitivity of the system to variations in the mtDNA degradation rate, where even a 2% change causes a huge difference in the mean and variance (Figure 4B), which is probably unlikely in vivo. This may hint at an underlying weakness of the model.
Reviewer #2:
This paper deals with a controversy in the mitochondrial inheritance field over the details of the inheritance bottleneck and the development of variability in mtDNA heteroplasmy levels across the cells of an organism. This new mathematical modeling approach is coming from a set of researchers, including one wellknown mitochondrial DNA expert, who were not involved in the original set of competing papers (primarily references Cree et al. 2008, Wai et al. 2008, Wonnapinij et al. 2010, Samuels et al. 2010, Cao et al. 2007, Cao et al. 2009). These authors also use simulation and mathematical methods that are far more sophisticated than the methods used in previous papers (references Cree et al. 2008 and Chinnery et al., 1999).
The paper does a good job of summing up the results and conflicting claims of the earlier literature. Their presentation of the differences in the assumptions of the competing models seems accurate to me, and they have built simulation models that encompass the three general competing sets of assumptions. The experimental data from the previous literature is all presented in a combined manner, which is valuable in itself. Through the Beyesian methods of assessing goodness of fit of the combined set of experimental data to the three competing models, the authors present a clear conclusion concerning the model that has the best agreement with the experimental data. This is a valuable contribution to the literature, and again it is very useful that this assessment is coming from an independent set of researchers.
While the mathematical detail that comes along with these more sophisticated methods is hard to digest, even for a trained mathematician, I would argue that it is essential to have this level of detail in the publication record.
Reviewer #3:
The authors have used stochastic modelling and Bayesian inference on published data sets to explore different models of the mitochondrial DNA genetics bottleneck.
There are several statements in the manuscript that are inaccurate and misleading. The notion that organisms ameliorate somatic mutation burden through bottlenecking is a hypothesis (not fact, as written), and not all oocytes have heteroplasmy as stated. Also, the idea that no model has been presented from the ‘bottom up’, with mtDNA as discrete elements, is not correct. The authors cite examples where this was the case (i.e. where the segregating units in the model are proposed to be mtDNA molecules). I also do not accept that there are no methods to allow the comparison of different bottlenecks. One of the authors (Poulton) is a coauthor on a Bayesian approach, and Samuels has published extensively on this (some of the papers cited in the current manuscript).
Some of the proposed mechanisms (which are modelled) have little or no experimental basis. For example, we do not know there is ‘turnover’ or mtDNA within germ cells, and certainly the rate is not known. We have no idea what the size of the nucleoid is in germ cells. There is no evidence that turnover is generally low during cell division to my knowledge, to name but a few. The problem is, by including all of these parameters, the potential outcomes of the modelling will increase.
Most importantly, their general conclusions are in keeping with the previous models developed by Samuels, which combine relaxed replication and random segregation contribute to the generation of heteroplasmy in the germ line (i.e. the Cree mechanism, as they describe it). In other words, their findings are confirmatory, not novel.
The theoretical manipulation of the bottleneck provides some interesting results, but these are somewhat obvious: if you reduce mtDNA content (by increasing degragadation; incidentally this is something not shown in the germ line), you will accelerate segregation and thus lead to the loss of mutations that reach high heteroplasmy levels. Likewise, it is not surprising that selective pressures influence heteroplasmy variance. Given that heteroplasmy levels are bound by 0% and 100%, any factors pushing the level in a particular direction will inevitably alter the range of values. It is also not surprising that the observed heteroplasmy variance is consistent with a range of bottleneck sizes, given the number of other parameters that can be varied in the described models.
https://doi.org/10.7554/eLife.07464.017Author response
Reviewer #1:
This is a resubmission of a pure modeling paper, with the addition of new experimental data which is presented as supporting the main conclusions of the modeling. The authors have carried out a very thorough theoretical and analytical analysis of the proposed mtDNA genetic bottleneck, and conclude that the existing stance taken by Cree and Wai model (b) is a sensible interpretation of the experimental data. As such, the current paper is reassuringly confirmatory. However, I still have concerns, in part relating to the way the paper is written, and in part relating to the new data, which raises additional questions.
Introduction. A statement is made: ‘No study yet exists combining an analytic ‘bottomup’ physical description of the bottleneck with mtDNAs as individual, discrete elements subject to a possible variety of replication and partitioning dynamics throughout an explicitly modelled series of cell divisions’. This is not correct, and I made this point in my earlier review of the paper. Their reference (Cree et al., 2008) includes such a ‘bottom up’ model. Although the compartments did not vary in that model, assumptions were made about the partitioning which are in accordance with the conclusions of this current paper.
Likewise, the statement ‘theory describing the behaviour of mtDNA populations throughout development, as opposed to estimating an overall effect of the bottleneck, is also absent’ is also an overstatement designed to justify the current work. There are several such theories that the authors refer to in their Introduction.
We appreciate the reviewer's guidance on the correct positioning of our paper. We do think that our work is a substantial and original contribution and agree that we must hit on the best way of conveying that contribution without inappropriately diminishing the excellent work of colleagues (we definitely do not want to do that!). The first statement in question contains several important qualifications which are not explicitly noted by the reviewer and which we suggest make our original wording acceptable:
“No study yet exists combining an analytic ‘bottomup’ physical description of the bottleneck with mtDNAs as individual, discrete elements subject to a possible variety of replication and partitioning dynamics throughout an explicitly modelled series of cell divisions” [our emphasis].
We would contend that this claim is in fact true. We are not claiming that ours is the first bottomup model of mtDNA dynamics; the novelty we claim is in the combination of our analytic progress (including features like explicit results for heteroplasmy variance with time) with the variety of possible replication and partitioning dynamics that we include. The simulations in Cree et al., 2008 do not provide analytic results nor consider other models than relaxed replication with binomial partitioning; previous analytic work has focused on particular specific models (for example, relaxed replication). We therefore contend that the novel combination of analysis and generality we claim are indeed true. We have reworded the text concerned to clarify that the novelty we claim lies in the variety of situations we consider (and thus the ability to select between these situations) and in the analytic progress we make.
We accept that our second sentence, regarding previous work modelling mtDNA through development, was still ambiguously phrased. We have simplified and just removed this sentence from the Introduction.
For the new experimental data, they use a line previously reported in their Cell Reports publication. A key message of this publication was that different mtDNA haplotypes behave differently, that selection occurs, and that this is influenced by the nuclear genetic background. All of these effects will influence the heteroplasmy variance. Although the authors touch on selection, they do not explore this thoroughly in the context of their previous work. How will this impact on their interpretation of their bottleneck models?
We thank the reviewer for making this important point. The possibility of selective pressures acting to change heteroplasmy in our models, in addition to the necessity of dealing with samples with different heteroplasmies, is the reason that we (and others) use the normalised heteroplasmy variance V'(h) throughout our work. This quantity corrects for the modulation of heteroplasmy variance V(h) by mean heteroplasmy E(h), and in so doing attempts to remove the effects of selective differences. This protocol has been used in previous work dealing with the NZB/Black mouse system which experience some (small) selective differences (Sharpley et al. Cell 151 333 (2012)). We have now drawn attention to this protocol in the main text. The influence of selective differences on observed heteroplasmy variance in our model can easily be computed by backtransforming the predicted normalised variance using the expected mean heteroplasmy, as can be seen from the structure of the transformation.
For the new experimental data, the heteroplasmy measurements are reported to a remarkable degree of accuracy. Do the authors really think that they can reliably report down to 0.01% as indicated in the supplementary data?
In reporting values to this precision we were following convention (e.g. Freyer et al. Nat Genet. 44 1282, 2012). The reviewer is correct in drawing attention to the fact that this reporting could be taken to falsely imply an accuracy of 0.01% in heteroplasmy levels. We have truncated the percentage figures given to 1 d.p. (a more reasonable representation of the data).
Previous theoretical work has shown that ∼50 heteroplasmy measurements are required to reliably compute heteroplasmy variance. Several of the published datasets do not achieve this, and thus the reported variance measurements are likely to be inaccurate. The same may also apply to the new experimental dataset.
Figure 4C raises concerns. For the ∼25 day old data, 6 data points (from a total of 15) are way outside the confidence intervals predicted by the various models. This implies that the new data set is not supporting the current or previous models.
The reviewer draws attention to the fact that sampling effects may mean that the variance of underlying heteroplasmy distributions is not adequately characterised by previous and current measurements. This is one reason we employ our particular inference protocol, which circumvents problems regarding the characterisation of a “true” heteroplasmy variance. ABC assigns support to a model based on its ability to recapitulate experimental observations, regardless of the interpretation of these observations. Each experimental observation is compared to a simulated outcome involving the same measurement protocol—so, for example, a heteroplasmy variance estimate of 0.01 derived from a sample of 20 cells is compared to 20 instances of stochastic simulation. Whether this estimate of 0.01 represents the “true” heteroplasmy variance does not matter for the inference process—and the distribution of the “true” heteroplasmy variance (corresponding to the variance of the underlying distribution from which samples are taken) is characterised by our approach.
The confidence intervals in Figure 4C do not correspond to the expected spread of individual datapoints, but rather to the spread of the mean variance. An appropriate analogy is to a traditional scatter of datapoints and s.e.m. bars, where many datapoints can lie outside the s.e.m. bars, because the spread of individual datapoints is broader than the spread of possible values for the mean. The fact that several datapoints lie outside the confidence interval for the mean does not therefore represent a weakness of our model. We have endeavoured to make this point clearer in the manuscript by including this discussion in the main text.
Finally, given the likely readership, it is important to keep reminding readers what is prediction, and what is actual measurement. The section on turnover illustrates this nicely. No turnover measurements have been made; indeed, they would be incredibly difficult to do.
Throughout our manuscript, experimentally determined data is plotted as circular datapoints and inferred behaviour as lines and shaded regions. We have made this clear in the main text at the point where we describe our first figure, and changed figure captions to explicitly include the fact that predictions have been made.
Reviewer #2:
There is clearly a great deal of variation in the new data for the normalized variance, but the new data does fit best with the authors' favoured birthdeathpartition model. The inclusion of this extra data, at the time points which are most discriminatory with regard to the various models, is therefore helpful, though not absolutely definitive. In particular, it looks as if the spread of the normalised variance itself has interesting dynamics, starting out low, increasing to a maximum at around 25 dpc before diminishing again. Is there any way in which the BDP model can capture these dynamics?
The apparently rich behaviour of variance is indeed intriguing, but due to technological limitations we cannot discount the possibility that this “variance variance” is due to sampling effects in individual datasets. Without technologicallyinaccessible large datasets we cannot confidently characterise higher moments like “variance variance” in any detail and instead focus on “mean variance”. Even this lowerorder moment, as discussed above, is subject to sampling issues, so satisfactory experimental characterisation of higher moments is currently impossible. The BDP model does not predict biphasic behaviour of higher moments; but the existing experimental data is not sufficient to characterise any departure between model and experiment as this level of detail.
Also it was not clear to me why the BDP curve in Figure 4A differs from the HB theory mean plot in Figure 4C.
The difference arises because the trace in Figure 4A arises from a single, bestfit set of model parameters, whereas the trace (and distribution) in Figure 4C describes the behaviour over the inferred posterior range on each parameter. For example, a given parameterisation may yield an onset of increasing V'(h) at a particular time (around 40dpc, in the case of the bestfit curve); but different parameterisations may have different values for this onset, “smearing” the increase over a range of times, as seen in the traces for averaged behaviour. We have described the difference between the traces in the main text.
Finally, I remain concerned that a 2% alteration in mtDNA degradation can have such a large effect, see Figure 5B. The authors must address this issue.
The reviewer raises an important rhetorical point which we have not yet made in a satisfactory way. In Figure 5A,B we display the effects of two different changes to the model parameterisation. In Figure 5A we change mtDNA degradation and apply a simultaneous and matching change in mtDNA replication. In Figure 5B we are attempting to display the effect of changing mtDNA degradation without a balancing change in mtDNA replication. Our previous attempt to display this effect involved perturbing the difference between degradation and replication rates by 2%. We then, perhaps confusingly, used this as a proxy for mitophagy.
We have, after consideration, used an alternative method. We now change the degradation rate itself by an additive constant, and keep the replication rate constant. This more directly represents the physical degree of freedom we want to address.
The reviewer comments that the magnitude of the effect provoked by a small change in degradation rate is large. This is indeed the case and is related to the distinction between Figure 5A,B above, and how mtDNA copy number is “controlled”. Copy number control in our simple model is achieved through replication and degradation rates that balance to produce the desired behaviour. The two rates can be varied in concert to maintain this behaviour while changing the timescale of mtDNA turnover, this is what we plot in Figure 5A. In Figure 5B we show the effects of perturbing this control by varying the rates independently. The effect of this “unbalanced” perturbation is then indeed pronounced, as we are applying a universal change to the system's expected behaviour. The comparable perturbation to a model applying control based, for example, upon a specific target copy number would be to vary the value of this target.
The selection of specific models for the control of mtDNA populations is a nuanced question and one that we are currently developing further mathematical approaches to address. Preliminary work in this topic has shown that for a range of different, plausibly parameterized, “control strategies” for mtDNA populations (including different incarnations of relaxed replication and setpoint controls, manifest through different birthdeath dynamics) the behaviour of heteroplasmy variance is very similar. We therefore believe that our claims regarding stochastic mtDNA turnover likely hold for a number of plausible scenarios, and that further work may shed light on the more specific details of cellular mtDNA control.
We have written more explicitly about these issues in the main text, underlining that our model employs this style of control and that other models (for example, those based on feedback control) could exhibit different behaviours. We will pursue the details of different control mechanisms in further work, but believe that our overall point—that increasing mitophagy may increase V(h) whether the system is tightly controlled or not—is valid and illustrated by our work.
[Editors’ note: the author responses to the previous round of peer review follow.]
Although there was enthusiasm for the quantitative approach to the mitochondrial DNA bottleneck, the reviewers felt that experiments directed at testing the predictions of the modeling are required for the manuscript to significantly advance the field. Therefore, at this stage, we are returning the manuscript to you with an encouragement to submit with experimental data. We hope the reviewers' comments are helpful.
We have taken your advice and tested predictions of our theoretical approach with an experimental collaboration. We now present an updated version of our study with novel experimental data from an mtDNA model system genetically different to those studied previously in this field. The new data confirms our theoretical predictions and further elucidates the mechanism of the mtDNA bottleneck, while also providing valuable insight about the generality of this process in genetically distinct systems. We believe that this inclusion indeed increases the power of our study, and hope that our improved manuscript will be considered for publication in eLife.
We have now, as suggested, collaborated with an experimental group to test the predictions of our theory. Encouragingly, these new experimental measurements demonstrate further support for the model we propose, in addition to demonstrating that these results are consistent across genetically diverse mtDNA pairings. The details of heteroplasmy distributions predicted by our model are also validated by this new data when we train our model on a subset of new datapoints and investigate the fit with the remaining, withheld subset.
On a broad level, Reviewers #1 and #2 were positive about our approach attempting to unify and amalgamate existing experimental data with consistent theory; Reviewer #1 and the editors suggested that confirmatory experimental results would strengthen our findings. We hope that the inclusion of new experimental data confirming our specific theoretical predictions addresses ensures the substance of our scientific contribution and also underlines the powerful cycle of using theoretical approaches both to amalgamate existing experimental data and to propose new experimental strategies to confirm theory and drive further scientific advance. Reviewer #3 raised the criticism that our selection of an existing theoretical model was not novel; we discuss the need for new models, as opposed to novelty in inference and in the analysis of existing models, below, but also hope that the new and confirmatory experimental data addresses some concerns of novelty.
Reviewer #1:
I do however have some concerns. Principally, the analysis is based entirely on previously published data, which is not all internally consistent. In particular, looking at Figure 2A for the variances, there doesn't seem to be very good agreement between the different datasets (especially comparing the Jenuth data with the rest). I appreciate that the authors' statistical analysis takes into account the experimental uncertainty, but the divergence between the different experimental data sets found in the different studies suggests that the effective experimental error bars may be much bigger than estimated. The ability to distinguish between the models via Bayesian model comparison may then be diminished. This point is important, as there is relatively little difference between the models with respect to the mean values (whose experimental values are also well clustered between the various data sets). Hence it does seem that the variances are the key to distinguishing between the models.
The reviewer highlights the discrepancies in the source data we use in our inference process. The measured values of heteroplasmy variance in these different studies do indeed cover a wide range. It is worth noting explicitly that these measurements are sample variances and, as the reviewer says, are subject to large uncertainties themselves. We are concerned with the behaviour of the underlying “expected variance”, which the measured datapoints can be envisaged as sampling. Depending on the sampling error involved (demonstrated under one model by Wonnapinij, Chinnery & Samuels (Am J Hum Genet 86:540 [2010])), all the existing data could be viewed as consistent with a given trajectory of this “expected variance”. Most measured variances will lie below this underlying value, as small samples typically do not capture the full variance of a population.
Our modelling approach captures this sampling error by recapitulating the corresponding experimental measurement protocol and comparing the measured outcomes. For example, if a given experimental datapoint reports a variance of 0.01 given 10 individual measurements, we draw 10 samples from the heteroplasmy distribution calculated from a given model and compare the variance of this sample against 0.01. By these means we align our computational simulation with the appropriate experimental protocol. Over the course of the inference and model selection processes, this process identifies those models that are likely to yield sampled variances most comparable to the data.
We have endeavoured to make this picture clearer in the text and in the caption to Figure 3. In addition, we have included new experimental data on heteroplasmy variance, obtained using a consistent experimental protocol that further and independently supports our results. We took these measurements at times chosen to best address the heterogeneous results in existing data and indeed find that the consensus picture from our theoretical treatment amalgamating previous, diverse data, is supported.
Also the statistical inference is based around a residual sum of squares difference. However, the way this was implemented (using logs of the copy number measurements, for example) seemed a bit arbitrary. Were other possibilities investigated to ensure consistency?
The reviewer notes that our chosen summary statistic involves some choices, including logs of copy number measurements. Our rationale here is to assign equal importance to matching the observed variability in copy number (which varies over orders of magnitude and is captured by many observations) and the observed variability in heteroplasmy variance (which does not vary over orders of magnitude and is captured by fewer measurements). The protocol employed weights the residuals arising from each of these classes equally, given the number of datapoints available for each. To address the reviewer's concern, we performed our model selection procedure with a range of other protocols, including varying this weighting over several orders of magnitude, and using logarithms of variance measurements as well as copy number. As now described in the manuscript, in all cases, our results held: the BDP model experienced the most statistical support.
More generally, I am concerned that the modelling is used primarily as a tool to fit existing data and thereby determine which model is most likely on the basis of previously published data. It would be more powerful if the authors could actually test the novel predictions of their modelling in new experiments. This would elevate the study to a higher level. Their model most definitely is predictive (Table 2), it's just that they ought to consider teaming up with an experimental group and prove these predictions right!
The reviewer expresses concern that our modelling approach is primarily a tool to fit existing data and perform model selection given available data, and encourages us to use it predictively. To address this point, we have taken the suggestion of the reviewer and the editors and performed experiments in a new mouse model of mtDNA heteroplasmy involving the pairing of a wildderived haplotype (labelled HB) with C57BL/6N.
These experiments show that, uninformed by the results in NZB/Balb, the bottleneck in a genetically different system is highly likely to manifest through the same mechanism we identify from previously published data. Moreover, the quantitative details of this mechanism are remarkably comparable between the two genetic systems, with similar amounts of variance increase due to copy number reduction and random turnover. We also test the predictions of our theory regarding the distributional features of heteroplasmy with time by training our model on a subset of these new data and find that our model successfully predicts the variance and distributions of the remaining data.
One other issue concerns the sensitivity of the system to variations in the mtDNA degradation rate, where even a 2% change causes a huge difference in the mean and variance (Figure 4B), which is probably unlikely in vivo. This may hint at an underlying weakness of the model.
The reviewer points out that our model is sensitive to the mtDNA degradation rate. The change in copy number and heteroplasmy variance provoked by a small change in degradation rate is indeed quite high. This theoretical perturbation is perhaps quite a “blunt tool”, corresponding to provoking a deliberate and systematic imbalance of mtDNA production and degradation. We may expect natural mechanisms to compensate for changes in degradation rate by invoking balancing changes in production rate. This is the picture in the previous subfigure, where turnover rather than degradation alone is altered, and the corresponding changes are less dramatic. We include the results for the “unbalanced” changing of mitophagy to aid intuition about how the model system behaves as independent degrees of freedom are varied. We do expect that more homeostatic, ‘meanreverting' mechanisms are likely to play a role in controlling mtDNA populations over longer timescales; however, the mathematical description of such models is more involved and so we work with this simpler case on the timescale of early development.
Reviewer #2:
The paper does a good job of summing up the results and conflicting claims of the earlier literature. Their presentation of the differences in the assumptions of the competing models seems accurate to me, and they have built simulation models that encompass the three general competing sets of assumptions. The experimental data from the previous literature is all presented in a combined manner, which is valuable in itself. Through the Beyesian methods of assessing goodness of fit of the combined set of experimental data to the three competing models, the authors present a clear conclusion concerning the model that has the best agreement with the experimental data. This is a valuable contribution to the literature, and again it is very useful that this assessment is coming from an independent set of researchers.
While the mathematical detail that comes along with these more sophisticated methods is hard to digest, even for a trained mathematician, I would argue that it is essential to have this level of detail in the publication record.
The reviewer noted the level of complexity of some of the mathematical analysis involved in the study. The full stochastic process treatment of the birthdeathpartitioning model is indeed rather complicated, though powerful; and indeed the derivation of common results used in the literature (such as expected heteroplasmy variance arising from Kimura's neutral model) is also complicated. However, the more straightforward results from this analysis, including expressions for mean heteroplasmy and heteroplasmy variance, are relatively straightforward to write down, and the full treatment affords a powerful way of reasoning about the dynamics of the system.
Reviewer #3:
There are several statements in the manuscript that are inaccurate and misleading. The notion that organisms ameliorate somatic mutation burden through bottlenecking is a hypothesis (not fact, as written), and not all oocytes have heteroplasmy as stated. Also, the idea that no model has been presented from the ‘bottom up’, with mtDNA as discrete elements, is not correct. The authors cite examples where this was the case (i.e. where the segregating units in the model are proposed to be mtDNA molecules). I also do not accept that there are no methods to allow the comparison of different bottlenecks. One of the authors (Poulton) is a coauthor on a Bayesian approach, and Samuels has published extensively on this (some of the papers cited in the current manuscript).
The reviewer highlights several misleading statements in the manuscript, specifically pointing out that:
i) The amelioration of mutational damage through the bottleneck is a hypothesis;
ii) Not all oocytes have heteroplasmy;
iii) Discrete models of mtDNA populations do exist;
iv) Statistical analyses of bottleneck mechanisms do exist.
It was not our intention to imply any certainty in points i and ii, and we have altered our Introduction and Abstract to be more compatible with these points. Our claim in point iii is not that no studies exist modelling mtDNA molecules as discrete elements (as the reviewer points out, several studies that we cite take this approach), but rather that no studies yet provide the combination of a broad and general bottomup model, analytic tractability, explicit timeseries modelling of mtDNA population size, and cell divisions. We have changed the appropriate text to make it clearer that we are referring to this combination. We have also altered the text to highlight the fact that statistical analyses of individual bottleneck models and individual datasets have been performed, and added a citation to the Bayesian study that the reviewer cites. We thank the reviewer for helping us clarify the manuscript.
Some of the proposed mechanisms (which are modelled) have little or no experimental basis. For example, we do not know there is ‘turnover’ or mtDNA within germ cells, and certainly the rate is not known. We have no idea what the size of the nucleoid is in germ cells. There is no evidence that turnover is generally low during cell division to my knowledge, to name but a few. The problem is, by including all of these parameters, the potential outcomes of the modelling will increase.
The reviewer states that several proposed mechanisms have little or no experimental basis, raising questions including the existence of mtDNA “turnover” and the size of the nucleoid in germ cells. They express concern that including these features may lead to a situation where “potential outcomes of the modelling will increase”, which we take to be a reference to the danger of overfitting—the fact that a suitably complex model can reproduce arbitrarily detailed behaviour but may not represent underlying “truth”.
As our work attempts to take a very general and allencompassing view of potential bottlenecking mechanisms, the inclusion of a set of features that have been postulated is important. Our Bayesian model selection approach automatically controls against overfitting; this is built into our sampling where the higher the number of parameters, the smaller is the probability that a perturbed set of parameter values will be accepted (analogous to the SMC case in Toni et al. J Roy Soc Interf 6 187 [2008]). As in standard Bayesian model selection, we are effectively integrating over all parameter values, and thus naturally penalising model complexity; unnecessarily complex models will be disfavoured. Furthermore, our parametric inference approach allows us to begin to address some of the uncertainties that the reviewer cites. Part of the message of the paper is that random turnover is the most likely explanation for measured heteroplasmy statistics in later development, and that a quantification of the magnitude of this turnover can be attempted, given data that currently exists. As another example, the reviewer states that we have no idea of nucleoid size in germ cells. But, as our results show, the inheritance of large (>∼ 5) clusters of mtDNAs has little support given existing experimental data, as such inheritance would provoke an increase in heteroplasmy variance rather higher than observed in existing measurements.
Finally, we of course accept (and explicitly state) that the model that we identify as the most supported does not represent unarguable truth. This is the motivation for the series of proposed experiments, by which the individual parameters that we consider can be further elucidated. Our proposed model does encouragingly match the new experimental results we include from a genetically distinct system, but more finegrained validation of the values of specific parameters should certainly be a target for future experimental work.
Most importantly, their general conclusions are in keeping with the previous models developed by Samuels, which combine relaxed replication and random segregation contribute to the generation of heteroplasmy in the germ line (i.e. the Cree mechanism, as they describe it). In other words, their findings are confirmatory, not novel.
The reviewer points out that our general conclusions are in keeping with models developed by Samuels, and argues that our findings are confirmatory rather than novel. An investigative study of existing hypotheses was indeed our intention with this paper, and we argue that unifying and amalgamating existing data to drive new theoretical progress does constitute a novel advance. The problem in the literature was not a shortage of models but an inability to identify the most relevant. We also point out that several of the insights that our work provides (see above) illustrate the power of a quantitative approach in making new theoretical statements, including restricting possible nucleoid sizes and showing the likelihood of mtDNA turnover.
Additionally, our mathematical approach is also a novel analytic solution using generating functions to fully and powerfully capture the stochastic behaviour of a natural system. Reviewers 1 and 2 point out that such a quantitative summary and treatment is valuable; we discuss the issue of novelty further above.
In response to comments from Reviewers 1 and 3 and the editors, we have now included novel experimental data from a genetically distinct system validating our theoretical results and supporting the quantitative predictions that our work makes. We hope that the combination of this new experimental data with the theoretical advances described above addresses concerns of novelty in our study.
The theoretical manipulation of the bottleneck provides some interesting results, but these are somewhat obvious: if you reduce mtDNA content (by increasing degragadation; incidentally this is something not shown in the germ line), you will accelerate segregation and thus lead to the loss of mutations that reach high heteroplasmy levels. Likewise, it is not surprising that selective pressures influence heteroplasmy variance. Given that heteroplasmy levels are bound by 0% and 100%, any factors pushing the level in a particular direction will inevitably alter the range of values. It is also not surprising that the observed heteroplasmy variance is consistent with a range of bottleneck sizes, given the number of other parameters that can be varied in the described models.
The reviewer claims that some of our results on manipulating the bottleneck are unsurprising. The qualitative results of some of the perturbations we investigate are indeed intuitive, including reducing mtDNA content to accelerate variance increase. We include these results for several reasons. First, as a quantitative characterisation of the statistically most supported model, illustrating the effect of varying different degrees of freedom of the model. Second, and more importantly, we include these results a quantification of these (perhaps qualitative intuitive) effects, making a numerical prediction regarding, for example, how much mtDNA reduction is required for a given increase in variance. Experimental tests of these predictions form part of our proposed strategies for further elucidation of the model.
https://doi.org/10.7554/eLife.07464.018Article and author information
Author details
Funding
Medical Research Council (MRC) (MR/J013617/1)
 Iain G Johnston
Biotechnology and Biological Sciences Research Council (BBSRC) (BB/D020190/1)
 Nick S Jones
Medical Research Council (MRC) (MR/J010448/1)
 Jo Poulton
Wellcome Trust (0948685/Z/10/Z)
 Jo Poulton
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Ethics
Animal experimentation: The study was discussed and approved by the institutional ethics committee in accordance with Good Scientific Practice (GSP) guidelines and national legislation. FELASA recommendations for the health monitoring of SPF mice were followed. Approved by the institutional ethics committee and the national authority according to Section 26 of the Law for Animal Experiments, Tierversuchsgesetz 2012  TVG 2012.
Reviewing Editor
 Jodi Nunnari, University of California, Davis, United States
Publication history
 Received: March 13, 2015
 Accepted: May 29, 2015
 Accepted Manuscript published: June 2, 2015 (version 1)
 Accepted Manuscript updated: June 4, 2015 (version 2)
 Version of Record published: July 1, 2015 (version 3)
Copyright
© 2015, Johnston et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,952
 Page views

 640
 Downloads

 63
 Citations
Article citation count generated by polling the highest count across the following sources: Scopus, Crossref, PubMed Central.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
Wound response programs are often activated during neoplastic growth in tumors. In both wound repair and tumor growth, cells respond to acute stress and balance the activation of multiple programs, including apoptosis, proliferation, and cell migration. Central to those responses are the activation of the JNK/MAPK and JAK/STAT signaling pathways. Yet, to what extent these signaling cascades interact at the cisregulatory level and how they orchestrate different regulatory and phenotypic responses is still unclear. Here, we aim to characterize the regulatory states that emerge and cooperate in the wound response, using the Drosophila melanogaster wing disc as a model system, and compare these with cancer cell states induced by ras^{V12}scrib^{/} in the eye disc. We used singlecell multiome profiling to derive enhancer gene regulatory networks (eGRNs) by integrating chromatin accessibility and gene expression signals. We identify a ‘proliferative’ eGRN, active in the majority of wounded cells and controlled by AP1 and STAT. In a smaller, but distinct population of wound cells, a ‘senescent’ eGRN is activated and driven by C/EBPlike transcription factors (Irbp18, Xrp1, Slow border, and Vrille) and Scalloped. These two eGRN signatures are found to be active in tumor cells at both gene expression and chromatin accessibility levels. Our singlecell multiome and eGRNs resource offers an indepth characterization of the senescence markers, together with a new perspective on the shared gene regulatory programs acting during wound response and oncogenesis.

 Cancer Biology
 Computational and Systems Biology
Colorectal cancer (CRC) remains a challenging and deadly disease with high tumor microenvironment (TME) heterogeneity. Using an integrative multiomics analysis and artificial intelligenceenabled spatial analysis of wholeslide images, we performed a comprehensive characterization of TME in colorectal cancer (CCCRC). CRC samples were classified into four CCCRC subtypes with distinct TME features, namely, C1 as the proliferative subtype with low immunogenicity; C2 as the immunosuppressed subtype with the terminally exhausted immune characteristics; C3 as the immuneexcluded subtype with the distinct upregulation of stromal components and a lack of T cell infiltration in the tumor core; and C4 as the immunomodulatory subtype with the remarkable upregulation of antitumor immune components. The four CCCRC subtypes had distinct histopathologic and molecular characteristics, therapeutic efficacy, and prognosis. We found that the C1 subtype may be suitable for chemotherapy and cetuximab, the C2 subtype may benefit from a combination of chemotherapy and bevacizumab, the C3 subtype has increased sensitivity to the WNT pathway inhibitor WIKI4, and the C4 subtype is a potential candidate for immune checkpoint blockade treatment. Importantly, we established a simple gene classifier for accurate identification of each CCCRC subtype. Collectively our integrative analysis ultimately established a holistic framework to thoroughly dissect the TME of CRC, and the CCCRC classification system with high biological interpretability may contribute to biomarker discovery and future clinical trial design.