Mitochondrial DNA heteroplasmy is modulated during oocyte development propagating mutation transmission

Description


Statistical tests
Statistical tests were performed on the transformed heteroplasmy data. The scope was to check whether there was a statistically significant difference in the heteroplasmy distributions over time. Transformed heteroplasmy is defined as [21] where h is the single cell value of heteroplasmy and h pup is an ear biopsy taken contemporaneously to the collection of oocytes. This definition of heteroplasmy has the advantage of showing the shift of the value h with respect to some initial value chosen (in this case, the aforementioned pup biopsy is the benchmark). Specifically, h positive (negative) corresponds to larger (smaller) values of h compared to h pup . Given how it is constructed, its sign and magnitude contain information about the sign and the strength of the selection. We designated all pre-antral cells as early and all antral cells late. We used a one-sample one-tailed t-test to explore shifts from zero: this allows us to reject the hypothesis that there is no-negative selection acting on the early cells. We used a two-sample one-tailed t-test to explore positive selection between the early and late cells. The tests assume a particular time as distinctive from the others: there are between 4 and 8 possible times we could have selected as distinctive. Both results remain significant after Bonferroni correction for 8 hypotheses.

Birth model
The dynamics of the system is simulated through a pure Poisson birth process encompassing the following reactions for mutants (m rs ) and wild-types (w rs ) subset of fixed size R 0 = m rs + w rs , whose composition (the proportion m rs to w rs ) varies in time. The ratio of m rs to w rs is defined to be the same as that in the whole cell: it is recomputed after each birth event. The idea of a fixed replicating subpopulation is implied in the immigration-like structure in [49] and yields a model with an approximately linear growth in population size (as is observed in the copy-number data). When cell heteroplasmy h exceeds a set threshold h th , the birth-rate of mutants is set to zero, while wildtypes continue to reproduce at a rate λ w . This is equivalent to introducing a mechanism of negative selection against mutants: this might be because, above the threshold, the replication is impeded or a process like selective mitophagy exists. If the total cellular mtDNA population is below R 0 we suppose that the dynamics is that of a conventional birth process.
In one of our models we suppose that a selective advantage only exists after a time T : before this time λ m = λ w , after this time they are free to take distinct values.
Also, in one of the models we consider that negative selection switches on at a value of the heteroplasmy that depends on the number of wild-types -specifically, at h th1 if w < w th , and at h th2 > h th1 if w > w th . This is consistent with Refs. including [16,51].
The parameters of the most general model are λ m , λ w , R 0 , h th1 , h th2 , w th , T (and H, a categorical parameter denoting which of six models is chosen, see below). For a given parametrization, we simulate the model using the Gillespie algorithm, also known as the stochastic simulation algorithm [50], where one stochastic trajectory corresponds to one cell in the experimental dataset.

Initial conditions
The experiment runs from day 1 to day 22 after the birth of the pup, so time will henceforth be measured in days post birth (dpb). Initial conditions for our model are the initial amount of mtDNA present in the cell n 0 and the initial heteroplasmy of the cell h 0 at day 1, where the heteroplasmy is defined as h = m/(m + w). Each measured cell can have a distinct initial conditionas inferred from the ear-biopsy of the corresponding pup.
The initial copy number is drawn from a Gaussian kernel density estimation using the data from the three distributions of the day 1 (corresponding to the three pups that contributed the oocytes). The distribution was truncated with a lower bound of n min = 75, which is the minimum value of the day 1 copy number data.
Our inference requires an initial heteroplasmy for each oocyte at day 1. We, instead, have knowledge of the ear-biopsy data for the pup from which the oocytes are removed. This ear-biopsy state is supposed to be the heteroplasmy value at an earlier time than day 1. We thus needed a model that takes as its input the pup ear biopsy for each pup and outputs a suitable candidate for its initial heteroplasmy at day 1. In order to build such a model, we assumed that the data that we already have from pups at day 1 contains the information about how selection acts on the initial (ear) heteroplasmy up to day 1. We thus first pooled together the day 1 transformed heteroplasmy distributions for each of the three pups obtaining h 0 . Then, we performed kernel density estimation (using a Gaussian kernel) to h 0 as shown in Fig.S4. The estimated density (violet in Fig. S4) serves as a baseline to obtain initial values for the transformed heteroplasmy. For pup i we can draw from the estimated density and build a vector h i ker with as many elements as the number of oocytes said pup contributed. Inverting (1) we can finally determine the initial heteroplasmy of oocyte j, pup i where h i pup is the value of the pup ear biopsy. In all instances the Gaussian kernel has been found the Scott's method to estimate kernel bandwidth.
An alternative approach is to simply suppose that the initial heteroplasmy of each oocyte at day 1 corresponds to the ear-biopsy heteroplasmy. Doing this in fact leads to very similar results.

Inference
Inference of the model parameters H, λ m , λ w , R 0 , T , h th,1 , h th,2 and w th was performed using the ABC rejection algorithm. H is an indicator parameter which encodes which model is being considered. We explored the following options: • H = 0, with positive and negative selection; • H = 1, with negative and without positive selection; • H = 2, with negative selection and positive selection starting at time T during development, whereas λ m = λ w before T ; • H = 3, without negative selection and positive selection starting at time T ; • H = 4, where neither positive nor negative selection are present; • H = 5, with positive selection throughout, but negative selection depends on the number of wild types: selection switches on at h > h th, Uniform priors are used for the heteroplasmy mitophagy thresholds h th and h th,2 ∼ U (0, 1), whereas the prior for h th1 is U (0, h th2 ). A log uniform distribution is chosen for the two replicative rates log 10 λ i ∼ U where E(h ) is the mean transformed heteroplasmy, E(n) is the mean copy number, V(h) is the heteroplasmy variance and the subscript specifies whether the statistics are calculated on the experimental (D) or simulated dataset (D ). We sum over all the pups except for the pups at t = 1 dpb which were used for setting the model's initial conditions. A j , j = 0, 1, 2, are normalization constants used to weight the three contributions approximately equally. Specifically, the summary statistics were normalized by the range of the corresponding experimental data. These are A −1 0 = 1.5, A −1 1 = 1.56 × 10 5 , A −1 2 = 1.47 × 10 −2 . We drew 6 × 10 6 samples from our prior for the ABC rejection algorithm.
The parametrizations were ranked according to their corresponding value of the distance metric ρ(D, D ), a smaller value of the latter corresponding to a stricter agreement to of the former to the data. We then looked at which model was more represented upon considering decreasing values of the distance metric, as a smaller threshold forces a stricter agreement with experimental data. As Fig. 2d (main text) shows, the model that best fits the data is H = 5, where negative selection depends on a threshold number of wild types. The best 200 parametrizations from model H = 5 have been retained as posterior distributions.
It is interesting to look at why the H = 5 model describes the data better than the others. If we look at the three error statistics we have considered (Fig. S5, data in red) we see that transformed heteroplasmy (panel (a)) dynamics seems to be split in two different phases: an early one and a later one, with a splitting point at some time in the range 14 − 20.  Fig. 3d (main text) we see that the worst performing models are those without negative selection (H = 3, H = 4). The poor ranking of these two models is likely due to there being very few data points with high heteroplasmy (h > 0.9): a feature that points to the existence of a purifying mechanism against very high levels of mutations (e. g. a negative selection like the one considered). To sum up, models 2 and 5 are both valid candidates in order to describe the trends in experimental data, but 5 performs better than the other one, and it is also more plausible from a biological standpoint [16,51].
In Figs. S6 and S7, we display approximate posterior distributions of the parameters from ABC for model H = 5. Small values for the size of the replicative subset R 0 (compared to population size of the cell) are favoured (Fig. S6c), because the exponential growth that the model undergoes until it reaches n = R 0 is incompatible with the linear copy number growth observed in the data. The posteriors of the birth rates λ m , λ w are better understood by looking at their covariation with R 0 (Fig. S8b,c) that clearly show the trade-off R 0 and their magnitude. This reflects the fact that the products R 0 λ m and R 0 λ w are peaked (mean ± standard deviation 6.3·10 3 ±1.1·10 3 ; 5.5· 10 3 ± 9.8 · 10 2 , Fig. S10a,b) because of the constraint due to copy number data that shows a clear linear trend (Fig. S5b). From the posteriors of the heteroplasmy thresholds h th1 , h th2 (Fig. S6d,e) (mean ± standard deviation 0.45 ± 0.20; 0.82 ± 0.1, respectively) it is clear that the negative selection is needed for the model to agree with the data (see discussion in sec. 3 for more details on this aspect). The distribution of the threshold wild-types value w th (Fig. S6f) at which negative selection changed intensity has the most support for values in the range 10 2 − 10 3 , which is early on in the dynamics. Lastly, since the relative advantage λ m /λ w > 1, the model fully supports mutants replicating faster than wild-types throughout the process (P = 99.5%), as expected.
To sum up, posteriors point to a model (H = 5) where below a moderatesized wild-type copy-number (w th = 10 2 − 10 3 ) the heteroplasmy threshold for tolerating mutants is low (h th1 ∼ 0.5) whereas above the characteristic wild-type copy-number the cells become more tolerant of mutant (h th2 ∼ 0.8).
The implication of this is that early in development, when copy-number is smaller, cells with higher initial heteroplasmy experience strong negative selection. This serves to pull down average heteroplasmy early in development. Later in development this negative selection is weaker (because the wild-types' copy-number is above w th ) allowing the continuous mutant positive selection to push up mutant-copy-number.  Figure S10: Posterior of the products R 0 λ w and R 0 λ m .

Analysis of a deterministic variant of the model
In order to better understand the dynamics of the chosen model (H = 5), we considered one of its deterministic variants. In this context the dynamics is described byṁ Equations (5) and (6) describe an exponential growth of mutants and wildtypes until a time t * when their total number equals R 0 , and a quasi-linear growth from t * onward. The two species replicate at their own rate λ m and λ w for mutants and wild types, respectively. Notice that the mutants' replication is blocked when h > h th,1 if w < w th , and when h > h th,2 if w > w th . From (5) and (6) it is easy to determine the rate of change of n and ḣ where r = λ m /λ w is the mutants' relative advantage.
The ODE system eqs. (7), (8) does not admit an analytical solution and was therefore simulated with different choices of the initial conditions and with parameters corresponding to the best parametrization 1 (Fig. S11). It is possible, though, to make some statements about the behavior of h and n after the exponential phase that is only a short transient and therefore does not determine the long-term dynamics. First of all we stress that the behavior of heteroplasmy and copy number depends on the relative value of initial heteroplasmy h 0 and the two thresholds h th,1 , h th,2 . Indeed, h 0 could be (a) larger than the two thresholds, (b) between h th,1 and h th,2 or (c) smaller than both 2 .
(a) h 0 > h th,2 > h th, 1 In this case heteroplasmy decreases and reaches h th,1 , then is constant until w = w th . After that, h increases again to hit h th,2 . Transformed heteroplasmy h is then always negative, but starts increasing after w = w th .
(b) h th,1 < h 0 < h th,2 Heteroplasmy behaves in quantitatively similar way to the previous scenario. The main difference is that h is negative when w < w th and increases to become positive afterwards.
(c) h 0 < h th,1 < h th,2 In this picture heteroplasmy increases to reach h th,1 and then remains constant until w < w th . After that, h keeps increasing again until it hits h th,1 , and again remains constant afterwards.
In all these cases n experiences a linear-like growth. Also, notice that since n ∼ t,ḣ ≤ λmR 0 n ∼ λmR 0 t , meaning heteroplasmy increases at a slower rate as the number of molecules grows. To understand why, it is sufficient to remember that the overall number of replicating molecules is fixed at a value R 0 , so while the total number of molecules grows over time, the ratio R 0 /n(t) decreases and so does the ability of the new molecules to affect the value of the heteroplasmy.
Let us now look at the numerical solutions of eqs. (7) and (8). In Fig. S11 we can see the time evolution of heteroplasmy, simulated in both the deterministic and the stochastic setting (in the latter case mean quantities are represented) with the parameters corresponding to the best parametrization, for different values of the initial heteroplasmy h 0 and with n 0 = 200. We can see that there is strong agreement between the two cases.