<figure><img src="image_1.png" /></figure>

HAL
open science

# The asymmetry of telomere replication contributes to replicative senescence heterogeneity

Thibault Bourgeron, Zhou Xu, Marie Doumic, Maria Teresa Teixeira

<figure><img src="image_2.png" /></figure>

## To cite this version:

Thibault Bourgeron, Zhou Xu, Marie Doumic, Maria Teresa Teixeira. The asymmetry of telomere replication contributes to replicative senescence heterogeneity. Scientific Reports, 2015, 5, pp.15326. 10.1038/srep15326. hal-01272075

HAL Id: hal-01272075

https://hal.sorbonne-universite.fr/hal-01272075

Submitted on 10 Feb 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_4.png" /></figure>

OPEN

# The asymmetry of telomere replication contributes to replicative senescence heterogeneity

Received: 22 June 2015

Accepted: 24 September 2015

Published: 15 October 2015

Thibault Bourgeron<fnref n="1" />,*, Zhou Xu²,*, Marie Doumic<fnref n="1" />,† & Maria Teresa Teixeira²,†

In eukaryotes, the absence of telomerase results in telomere shortening, eventually leading to replicative senescence, an arrested state that prevents further cell divisions. While replicative senescence is mainly controlled by telomere length, the heterogeneity of its onset is not well understood. This study proposes a mathematical model based on the molecular mechanisms of telomere replication and shortening to decipher the causes of this heterogeneity. Using simulations fitted on experimental data obtained from individual lineages of senescent *Saccharomyces cerevisiae* cells, we decompose the sources of senescence heterogeneity into interclonal and intraclonal components, and show that the latter is based on the asymmetry of the telomere replication mechanism. We also evidence telomere rank-switching events with distinct frequencies in short-lived versus long-lived lineages, revealing that telomere shortening dynamics display important variations. Thus, the intrinsic heterogeneity of replicative senescence and its consequences find their roots in the asymmetric structure of telomeres.

Heterogeneity in cellular processes, whether programmed or stochastic (also called noise), has important physiological consequences in all organisms. In bacteria and microorganisms in general, phenotypic variability generally increases the fitness of the population both in normal and adverse conditions<fnref n="1" />⁻⁴. Developmental and, more generally, differentiation programs also aim at producing heterogeneity out of a uniform population of cells or even a single cell⁵⁻¹⁰.

In addition to the diversity of states reached by differentiation, many terminal differentiation processes also show important variations in time, from the first differentiation signal to the final cell fate commitment⁵,¹¹⁻¹³. In unicellular organisms, this temporal heterogeneity may be viewed as a survival strategy by promoting flexibility in fluctuating environments. The fundamental role of noise in all these processes is further confirmed by the fact that noise can be controlled and evolved in genetic circuits to optimize its functionality¹⁴⁻¹⁶.

Replicative senescence is an asynchronous process well documented to display interclonal and intraclonal variations in different organisms¹⁷⁻²⁰ and first described as a limit to the proliferation of human primary fibroblasts, the Hayflick limit²¹. An important step in understanding senescence onset and heterogeneity was to relate the Hayflick limit to the molecular mechanism of telomere shortening. Indeed, the ends of eukaryotic linear chromosomes, called telomeres, require a specialized mechanism for its maintenance, in most cases the holoenzyme telomerase. Without maintenance, as in primary fibroblasts

<footnote n="1">¹INRIA Paris-Rocquencourt, Domaine de Voluceau, BP 105, 78153 Le Chesnay, France. UPMC University of Paris 6, JL Lions Lab., 4 place Jussieu, 75005 Paris, France. ²Sorbonne Universités, UPMC Univ Paris 06, CNRS, UMR8226, Laboratoire de Biologie Moléculaire et Cellulaire des Eucaryotes, Institut de Biologie Physico-Chimique, 75005 Paris, France. ³These authors contributed equally to this work. ⁴These authors jointly supervised this work. Correspondence and requests for materials should be addressed to Z.X. (email: zhou.xu@ibpc.fr) or M.D. (email: marie.doumic-jauffret@inria.fr)</footnote>
or in telomerase-negative *Saccharomyces cerevisiae* cells, telomeres shorten and induce replicative senescence after many cell divisions. Many other factors affect and modulate replicative senescence and its inherent variability, including oxidative stress<sup>22-24</sup>, mitochondrial dysfunction<sup>25,26</sup> and environmental stress<sup>27</sup>, which may act together and/or through telomeres<sup>23,28,29</sup>. Senescence may thus be seen as the outcome of complex interactions between metabolic and DNA-damage (including telomere dysfunction and attrition) pathways<sup>30</sup>, contributing to increasingly more variability over time or cell divisions<sup>31</sup>. This stochastic view of senescence stands in contrast with the idea that telomere shortening may represent a molecular clock mechanism for cell divisions. It is thus critical to assess the contribution of telomere-related variations to senescence heterogeneity.

In cells that express telomerase such as stem cells or most cancer cells in mammals, or unicellular eukaryotes, telomere length is a dynamic equilibrium between shortening and elongation by telomerase. Thus, in the presence of telomerase, the telomeres in a single cell encompass a broad range of lengths and an individual telomere, even in a clonal population of cells, displays an important length heterogeneity<sup>32</sup>. Many studies in different organisms converge to the idea that, when telomerase is absent, the average length of all telomeres is not critical for senescence onset. Rather, the shortest telomere(s) reaching a threshold may be the actual signal for senescence arrest<sup>33,34</sup>. Therefore, the length heterogeneity of the shortest telomeres, which may be greater than that of the others<sup>20</sup>, should strongly impact on senescence onset heterogeneity. At the molecular level, telomeres end with a 3'-end single-stranded DNA overhang<sup>35,36</sup>. Such a structure stems from the inability of the replication machinery to integrally replicate telomeres<sup>37,38</sup>. The direct consequence of this end-replication problem is the asymmetry of telomere replication, resulting in their overall shortening when no maintenance mechanism is present<sup>39</sup>. But whether this shortening mechanism by itself generates heterogeneity in telomere length has not been specifically addressed.

Here we aim at characterizing the different sources of senescence heterogeneity by a computational approach, using *Saccharomyces cerevisiae* as a model organism for telomere-induced senescence. In a recent study, we reported and quantified replicative senescence heterogeneity at the level of individual lineages, which was a first step to characterize heterogeneity by overcoming average measurements in population studies<sup>40</sup>. To analyze the heterogeneity observed in that study, we model the asymmetric replication of each telomere in a cell, using recent data on the molecular mechanism of the 3'-end overhang processing<sup>41</sup>, and the random segregation of the chromosomes in the two daughter cells. We consider senescence signal as the result of one or more telomeres reaching a length threshold and thus focus on telomere-dependent mechanisms of generating heterogeneity. By comparing simulations based on the mathematical model with the experimental data, we could decompose senescence heterogeneity into fundamentally distinct parts at different levels and evaluate their contribution.

We find that, in the model system, not only do cell-to-cell variations in the initial telomere length distribution contribute to senescence heterogeneity, but the molecular mechanism of asymmetric telomere replication also generates important intrinsic variations in telomere length dynamics. Therefore, a substantial part of replicative senescence heterogeneity may be structurally built in the telomere replication mechanism. As telomere structure and replication mechanism are broadly conserved in eukaryotes with some notable exceptions, this study provides a new molecular basis for phenotypic variations in multicellular organisms, with pleiotropic consequences in differentiation, aging and cancer emergence.

## Results

**A model of telomere replication and its impact on senescence onset.** We define a lineage as the consecutive cells with a mother/daughter relationship, which we were able to track in a dedicated microfluidics-based live-microscopy assay<sup>40</sup>. The aim of the following model is to track the 32 telomeres of a haploid yeast cell at the single-cell level in the consecutive cells that make up a lineage, in the absence of telomerase. We note $L = (L_1^n, ..., L_{32}^n)$ the lengths of the 32 telomeres of a cell at the n-th generation in a lineage.

Telomere length dynamics is different in the presence and in the absence of telomerase. We first describe the steady state telomere length distribution obtained in the presence of telomerase, which we use as the initial distribution for all subsequent simulations (Fig. S1). We then describe a model based on the molecular mechanism of telomere replication and shortening in the absence of telomerase in individual lineages.

*Initial steady state telomere distribution in the presence of telomerase.* In the presence of telomerase, the telomeres in a cell evolve by two means: (i) they are shortened on average after replication, and (ii) they can be elongated by telomerase with a probability that decreases with increased length<sup>42</sup>. A model for the evolution of each length $L_k^n$, for $k = 1, ..., 32$, was studied in<sup>20</sup>:

$$L_k^{n+1} = \begin{cases} L_k^n - s, & \text{with probability } 1 - P(L_k^n) \\ L_k^n - s + B, & \text{with probability } P(L_k^n), \end{cases} \quad (1)$$

where s is the mean shortening rate at each division, B is the random elongation length due to the action of the telomerase and P(L) is the probability of action of telomerase as a function of telomere length. The
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Model of telomere replication and segregation. (a) The molecular mechanism of telomere replication, at the basis of the mathematical model, is depicted. Replication of the parental telomere, harboring a 5–10-nucleotide 3'-end overhang, leads to two new telomeres of different lengths through the leading and lagging strand synthesis: one with parental length (black and green) and the other with shortened length (shorter by the overhang length; grey and red). Random segregation of chromosomes into the two daughter cells shuffles the set of 64 telomeres with each daughter inheriting 16 telomeres of parental length and 16 shortened ones. (b) Scheme of a branching model for cell division. Only one of two daughter cells at each generation is kept in the simulations as in the experiments, thus defining a lineage. Telomere length distribution is recorded for each cell in a lineage.</figcaption></figure>

variables s, B, P(L) do not depend on k or n. As in$^{20}$, the parameter s is set at s=3.5 bp, B is a random variable following a geometrical law of parameter p = 0.026, and lastly the probability P(L) is fitted from the cumulative distribution function of empirical elongation lengths$^{42}$, and expressed as follows:

$$P(L) = \begin{cases} (1 + \beta(L - L_0))^{-1}, & \text{if } L \ge L_0 \\ 1, & \text{if } L < L_0. \end{cases} \qquad (2)$$

This probability of action of telomerase depends on two parameters β and L₀ that can change from a yeast strain to another and a condition to another, depending on the *in vivo* activity of telomerase. We assume that the telomerase-positive initial cells we randomly pick for simulations have undergone a large number of divisions, thus reaching a steady state for their telomere lengths. Under this assumption, we can fit β and L₀ from only two parameters of telomere length distribution: the mode and the skewness. The mode of telomere length distribution may be experimentally measured by telomere restriction fragment (TRF) Southern blot and the skewness is assumed to be the same as in$^{20}$. For more details, see Materials and Methods. We are thus able to simulate steady state telomere lengths in a cell, corresponding to the experimental strain TetO2-TLC1 yT528 used in$^{40}$, which had a mode of 260 ± 16 bp. This steady state distribution is used in our simulations as the starting point distribution before telomerase inactivation (Fig. S1).

**Telomere replication and shortening in the absence of telomerase.** In the absence of telomerase, the average telomere length decreases at each division. At the level of individual telomeres, the mechanism of shortening stems from the DNA-end-replication problem$^{37–39,41}$. In late S-phase, DNA replication gives rise to 64 telomeres: 32 lagging-strand telomeres of parental length and 32 leading-strand telomeres with shortened lengths$^{41}$ (Fig. 1a). After random segregation of the two sets of chromosomes, each daughter cell receives 16 telomeres of parental length and 16 shortened telomeres, because one extremity of a chromosome is the result of the leading strand synthesis while the other is replicated by the lagging strand machinery. For the shortened telomeres, the loss of nucleotides corresponds to the length of the overhang of the original telomeres before replication, which is 5–10 nucleotides$^{41}$ (Fig. 1a). We thus write the following mathematical model to describe telomere replication and segregation into daughter cells: for a random permutation σ of {1, ..., 32}, for k = 1, ..., 32,

$$L_k^{n+1} = \begin{cases} L_k^n, & \text{if } k \in E \\ L_k^n - s_k, & \text{if } k \in E^c, \end{cases} \qquad (3)$$

where E = σ({1, ..., 16}), E^c = σ({17, ..., 32}) and s_k are uniform random variables on the set {5, 6, 7, 8, 9, 10}, corresponding to the overhang length. At generation n + 1, we keep track of only one of the two daughters, the succession of cell divisions thus defining a lineage (Fig. 1b). Another hypothesis was tested for the law of the random variables s_k, and will be explained below. For n ≥ 0, the smallest
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Fit of senescence onset timings. (a) The type A lineages from$^{40}$ were used as experimental data. Display of n=24 independent TetO2-TLC1 cell lineages. Each horizontal line represents a lineage, and each segment is a cell cycle. Cell-cycle duration is indicated by the colour bar. The timings of transition into a senescent state were extracted (black line). (b) Representative set of 5 simulations of 24 senescence timings in independent lineages (blue lines, Lmin = 19 bp, a = 0) compared to the experimental senescence onset timings (black line), extracted from (a). (c) Envelope of n = 1000 simulations (blue-shaded area). (d) 3D-plot of the error of the fit *e* as a function of *a* and *Lmin*. (e) 2D-plot of the error of the fit *e* as a function of *Lmin* for *a* = 0.</figcaption></figure>

value among $L_1^n, ..., L_{32}^n$ is denoted by $L1$ and the second smallest value is denoted by $L2$. We model senescence onset as the consequence of a signal triggered by the weighted contribution of the two shortest telomeres reaching a threshold $Lmin$. Thus, the cell that stops dividing verifies the condition $L1 + aL2 < Lmin$ for a constant value *a* in [0, 1]. The generation of senescence is then recorded. The dynamics of the shortening of the 32 telomeres in the successive cells of lineage is illustrated in the Supplementary Movie.

**The shortest telomere is the sole contributor to senescence onset signal.** As the model was built to analyze telomere length dynamics at the single-cell level, we used equally high-resolution data, namely senescence onset timing of individual yeast cell lineages, to validate it$^{40}$. In that study, we clus- tered telomerase-negative lineages into two types depending on the presence of precocious and reversible cell cycle arrests before transition into senescence. The majority of the lineages, called type A, did not display these arrests and instead underwent normal cell divisions until a single sharp transition into senescence (Fig. 2a). Their progression into senescence was consistent with a model of progressive tel- omere shortening and sudden signal for senescence arrest. We thus extracted the experimental profile of senescent onset timings from the type A lineages and used these data to fit our model (Fig. 2a-c).

Evidence suggest that a single critically short telomere is sufficient to signal senescence$^{34}$. In the model, the expression of senescence trigger as $L1 + aL2 < Lmin$ was chosen to test this hypothesis by evaluating the contribution of the second shortest telomere to senescence signaling. We ran a total of >26 million simulations with different values for *a* and *Lmin* and found that the parameters that best fitted the experimental results were *a* = 0 and *Lmin* = 19 bp, with a root mean square deviation, which we simply call error of the fit, of *e* = 0.59 (Fig. 2d,e). An example of a fit with an error *e* = 1.7 (*a* = 0 and *Lmin* = 60 bp) is shown for reference in Fig. S2. This indicated that a contribution of the second shortest telomere to senescence signaling was unlikely and the shortest telomere reaching a length threshold was a sufficient condition to fit the data.

To test the robustness of the fit and assess whether the fit errors for *a* > 0 were significantly different from the minimal error found for *a* = 0, we changed the initial mode of telomere length distribution and tested $260 + 16 = 276$ and $260 - 16 = 244$ bp, corresponding to the upper and lower limits of the uncer- tainties of the initial telomere length (Fig. S3)$^{40}$. Using these limit values, we found *a* = 0 for both condi- tions and *Lmin* = 35 and 2 bp, respectively (Fig. S3). Importantly, the errors in both cases (*e* = 0.60 and
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Telomere shortening dynamics at the population level. (a) Representative TRF Southern blot showing telomere shortening after telomerase inactivation. Xho I-digested terminal restriction fragments (TRF) were probed with a radioactive telomeric probe and average telomere length was measured.</figcaption></figure>

(b) Comparison between experimental measurements of the average telomere length as in (a) ($n=3$ independent experiments, the error bars represent the standard deviation) with simulation of the dynamics of the average telomere length in a population of telomerase-negative cells assuming $Lmin = 19$ bp. The upper and lower limits of the simulation defining the uncertainty region (grey-shaded area) were derived from the cases $Lmin = 35$ and $2$ bp, respectively.

$e=0.60$) were the same as for the first fit ($e=0.59$), discarding the possibility that the higher values of the error found for $a > 0$ could be non significant considering the uncertainties limits of the initial mode of telomere length distribution experimentally measured ($260 \pm 16$ bp) and used as an input parameter.

We also estimated the lower bound of the error *emin* and the typical error *etyp* between one simulation and all the others. The variance of $N=1000$ simulations immediately yielded $emin = 0.22 \pm 0.005$, which represents the intrinsic variability of the simulated data. However, *emin* does not provide an estimate of the expected error between a specific set of data (whether experimental or simulated) and the bulk simulations. We thus calculated the error of one simulation compared to the others, which is formally the same calculation as for *e* but with simulated data instead of experimental ones, and we found *etyp* = $0.35 \pm 0.10$. We notice that the errors of the fits with our model are significantly greater than *etyp* suggesting that there may be some mechanisms that we have not taken into account in the model. Nonetheless, because the experimental data was completely enclosed within the envelope of our simulations (Fig. 2c), the model still captured the major mechanisms driving replicative senescence.

To take the precision of 3'-end overhang measurement into account⁴¹, we also tested a shorter overhang length (5-7 nucleotides) in the model and found $a=0$, $Lmin=47$ and $e=0.45$ (Fig. S4), suggesting that a lower shortening rate of 3 bp/division might be slightly more consistent with the experimental data. However, even here, no contribution of the second shortest telomere could be found.

Taking all these fits together, the values of *a* and *e* were robust, in contrast to *Lmin*, which fluctuated depending on the input parameters. Thus, our model robustly reproduced the experimental timings of senescence onset, notably its heterogeneity, and reinforced the idea that only the shortest telomere contributed to signal senescence.

**Average telomere shortening in a population of senescent cells.** In a population of senescent yeast cells, the average telomere length would reach over time a plateau at 80–140 bp¹⁸ (Fig. 3a). Given this average value for all telomeres, a threshold for the shortest telomere at $Lmin = 19$ bp (or at least between 2 and 35 bp) may not be intuitive because it is not trivial to derive the length of the shortest telomere from the large spread of telomere length distribution. To challenge our model, we predicted the shortening dynamics of the average of all telomeres by simulations with our model using a single threshold constraint on the shortest telomere and compared it to experimental measurements.

To do so, we simulated the divisions of a population of cells using a simple branching model, starting with a mean telomere length of 260 bp. As type A cells display homogeneous and constant cell cycle durations until senescence onset (Fig. 2a)⁴⁰, we chose a binary behavior for simulated cells, being either in a dividing or an arrested state. A cell arrests when its shortest telomere reaches the $Lmin = 19$ bp threshold.
At different time points, the length of all the telomeres in the population was averaged and recorded. This procedure mimicked the experimental assay in which telomere length was measured by TRF Southern blot. The simulated average telomere length decreased linearly over the first ~40 population doublings with a rate of ~3.5bp/generation and then stabilized to reach a plateau at ~130bp (Fig. 3b). This result closely matched the experimental telomere length dynamics using a telomerase-inactivated strain and TRF analysis (Fig. 3b). The threshold value of Lmin = 19bp used for the shortest telomere in this simulation was thus consistent with and validated by biological data.

Strikingly, in the simulation, the threshold rule for the shortest telomere alone could constrain the dynamics of the average telomere length of the whole population of cells. Indeed, the average telomere length reached a plateau because of a competition effect: cells that reached senescence later were constantly selected over time and would have an average telomere length in the range of the plateau value. This competition effect reflects the complex dynamics of telomere shortening leading to a heterogeneous mixture of cells as they progress toward senescence.

**Interclonal variations partially contribute to senescence heterogeneity.** We formulate two hypotheses to explain the experimental heterogeneity. First, heterogeneity in senescence timing may be explained by interclonal variations of the shortest telomere length in the founder cells, in which telomerase was inactivated. The other hypothesis involves an intrinsic variability of telomere length generated over the course of senescence. To discriminate between these two possibilities, we reasoned that the heterogeneity due to cell-to-cell (interclonal) variations would be suppressed if a single founder cell, instead of different independent ones, were used to produce multiple senescent cell lineages. In contrast, the intrinsic variability of senescence onset would be the same in both cases.

We thus simulated multiple cell lineages using the same initial cell, with its specific telomere length distribution, and analyzed the heterogeneity of the subsequent senescence onset profile by calculating the coefficient of variation (CV), defined as the standard deviation of the set of senescence onset timings divided by their mean (Fig. 4a,b). As expected, the simulation with a single founder cell did not fit the experimental data well because the actual experiment used different independent founder cells (Fig. 4d, compare with Fig. 4c). We calculated the coefficient of variation and found $CV = 0.16 \pm 0.05$. In striking contrast, the simulations in Fig. 2, in which 24 independent founder cells were used, corresponded to a greater coefficient of variation of $0.24 \pm 0.04$ respectively (Fig. 4c). The decrease of CV upon using a single founder cell in the simulation showed that interclonal variations of telomere lengths contributed to the heterogeneity of senescence onset. However, we could still appreciate substantial variations in senescence onset in the clonal experiment simulations, suggesting the existence of yet another source of intraclonal noise (Fig. 4d), which seemed to be intrinsically captured in our model.

**The asymmetry of telomere replication generates stochasticity in telomere shortening dynamics.** The only mechanism implemented in our model that might lead to differences in daughter cells regarding senescence onset would be the telomere strand asymmetry due to the 3'-end single-stranded DNA overhang and the distinct machineries involved in the replication in the two strands, namely the leading and lagging strand replication machineries (Fig. 1a). Indeed, asymmetric replication of telomeres produces two new telomeres of different lengths (one equal to the parental telomere, the "lagging telomere"; the other shorter by 5-10bp, the "leading telomere";⁴¹) so that, at the level of an individual telomere and cell division, telomere shortening is not constant but probabilistic with an all or nothing pattern.

To investigate the contribution of asymmetric telomere replication to senescence heterogeneity, we slightly modified our model by forcing, for each telomere replication, the two newly replicated telomeres to have the same length, which is the parental length minus the average shortening rate of the standard model. We then tested this artificial symmetric telomere replication model starting with a single founder cell (Fig. 4f). We found virtually no variations in senescence onset for all the lineages ($CV = 0.04 \pm 0.01$, compared to $CV = 0.16 \pm 0.05$ in the case of asymmetric telomere replication), showing that we managed to remove nearly all sources of variations. The decrease in the coefficient of variation with this symmetric model could barely be observed when we simulated independent founder cells, as it was largely masked by the interclonal variations ($CV = 0.20 \pm 0.03$, compared to $CV = 0.24 \pm 0.04$ for the asymmetric model) (Fig. 4e). Thus, in our mathematical model, the asymmetry of telomere replication directly affects senescence heterogeneity by modulating telomere-shortening dynamics in distinct lineages, even starting with the same founder cell.

**Initial identity of the shortest telomere at senescence onset.** Theoretically, by this asymmetric replication, the shortest telomere is maintained in one of the two daughter cells and shortened in the other. In consequence, in some rare lineages, the shortening rate of the shortest telomere may be much lower than average and in others, much higher, thus explaining the intrinsic variations. However, the picture may be more complex: the shortening dynamics of the other telomeres has to be taken into account. For example, in a specific lineage, if the lagging telomere of the initial shortest telomere is inherited for several divisions in a row, the second shortest telomere may reach at some point the same length as the shortest and eventually switch rank with it and become the "new" shortest (See examples in Supplementary Movie). Indeed, stating that the shortest telomere signals senescence does not necessarily
<figure><img src="image_8.png" /><figcaption>c Asymmetric telomere replication model</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>e Symmetric telomere replication model</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>d Asymmetric telomere replication model</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>f Symmetric telomere replication model</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 4.** Contribution of interclonal variations and the asymmetry of telomere replication to senescence heterogeneity. (a) The scheme for the simulations used in (c,e) is pictured, showing independent founder cells for each lineage. (b) The scheme for the simulations used in (d,f) is shown and depicts a single founder cell from which multiple lineages are derived. (c) Representative set of 5 simulations starting from independent founder cells and using a model of asymmetric replication of telomeres (standard model). For (c–f), the mean and standard deviation of the coefficients of variations are indicated (CV, *n* = 55 independent simulations of 24 lineages). (d) Representative set of 5 simulations starting from a single founder cell and using a model of asymmetric replication of telomeres. (e) Representative set of 5 simulations starting from independent founder cells and using a model of symmetric replication of telomeres (modified model). (f) Representative set of 5 simulations starting from a single founder cell and using a model of symmetric replication of telomeres.</figcaption></figure>

imply that this specific telomere was *initially* the shortest telomere when, many divisions earlier, telomerase was inactivated.

To quantify the frequency of these “rank-switching” events, we used the simulation shown in Fig. 2 and traced back the origin of the signaling shortest telomere to find its initial rank in the founder cell, in which telomerase had been inactivated. Unexpectedly, in nearly 40% of the lineages, the initial shortest telomere was not the shortest at senescence (Fig. 5a). For instance, in ~20% of the cases, the initial second shortest telomere eventually became the shortest (Fig. 5a).

We then investigated the dependence of this telomere rank-switching phenomenon on the timing of senescence onset (Fig. 5b). We found that the short-lived lineages were essentially determined by their initial shortest telomere whereas for longer lineages, the initial second, third, etc., contributed increasingly more to signal arrest at the transition to senescence (Fig. 5b, Supplementary Movie). Nevertheless, the initial shortest telomere always maintained a major role in signaling senescence (Fig. 5a,b), consistently with previous work²⁰. This simulation result suggested that *in vivo* short-lived and long-lived senescent lineages may experience different telomere shortening dynamics, again highlighting the complexity and heterogeneity of senescence dynamics.
<figure><img src="image_14.png" /><figcaption>**Figure 5.** Telomere rank switching events involving the shortest telomere. (a) Aggregate frequency of rank switching events, i.e. frequency at which the initial $n^{th}$ shortest telomere eventually becomes the shortest at senescence onset. Data derived from 10 000 simulations. (b) Frequency of rank switching events as a function of the generation of senescence onset. The initial telomere rank is indicated by its color.</figcaption></figure>

Discussion

Replicative senescence is a highly asynchronous and heterogeneous process. Characterizing the source of the underlying heterogeneity has important and broad implications that range from fundamental understanding of senescence signal and the escape from it in survivors or cancers to improvement of diagnosis and treatment of telomere-related diseases. Previous mathematical models of telomere shortening focused on proliferation potential and on factors that regulate telomere state and length, in populations of cells and mostly in mammalian cell cultures, in which complex pathways interact to control cell proliferation<sup>22,29,43-52</sup>. In this study, we break the overall senescence heterogeneity down into mechanistically distinct components in an arguably simpler model for telomere-induced replicative senescence, *S. cerevisiae*. To do so, we modeled the telomere shortening mechanism at the molecular level using the most accurate measurement of the 3'-end overhang, taking the recent experimental demonstration of the asymmetry of telomere replication into account<sup>41</sup>, and its consequences on senescence signaling at the level of individual cell lineages. Simulating and analyzing single-cell lineages instead of populations of cells were necessary to accurately assess senescence heterogeneity and we took advantage of single-lineage data of telomerase-negative yeast cells undergoing replicative senescence to fit and support our simulations<sup>40</sup>.

By comparing the simulations based on the model to the experimental data, we could test hypotheses on the nature of senescence signal. We found that there was robustly no contribution of the second shortest telomere to senescence, quantitatively reinforcing the idea that the shortest telomere reaching a threshold is sufficient to signal senescence<sup>20,34</sup>. We could also show by simulations that this rule constrains the division dynamics of senescent cells in a population, thus explaining the two-phase shortening pattern of the average telomere length dynamics<sup>18</sup> (Fig. 3b).

After establishing the nature of senescence signal, we went on characterizing the different sources of heterogeneity. Clearly, cell-to-cell variations in the initial telomere distribution, due to the dynamics of shortening and lengthening by telomerase occurring independently at each telomere<sup>32</sup>, contribute to senescence heterogeneity but only partially. These interclonal variations were suppressed when we analyzed by simulations the progeny of a single founder cell instead of lineages from independent founder cells. Strikingly, in this condition, we still observed important variations in senescence onset timings and we could attribute them to the asymmetry of telomere replication mechanism. Indeed, this asymmetry is a molecular source of noise leading to probabilistic telomere shortening at the level of individual telomeres. Consequently, individual telomere shortening rate may be different in different cell lineages. Furthermore, not only does this mechanism modulate the shortening rate of the initial shortest telomere, but it also leads to telomere rank-switching events (Supplementary Movie). To illustrate how strongly this mechanism affects senescence dynamics, let us consider two opposite examples. In the specific cell lineage that systematically inherits the strand replicated by the leading-strand machinery of shortest telomere, thus shortening by ~5-10 bp at each division, replicative senescence may occur as fast as ~10 divisions. In contrast, in the long-lived lineages, the successive cells often inherit the G-strand of the shortest telomere, thus maintaining its length. But eventually, the second shortest (and afterwards even the third, fourth, etc.) reaches the length of the shortest and switches rank with it. Therefore, in the long-lived lineages, the two or three shortest telomeres tend to shorten at a slow rate, and the cell lineage can divide for more than 50 divisions until senescence. Overall, short-lived and long-lived lineages may experience very different telomere shortening dynamics and thus senescence dynamics, with potential consequences in the mechanisms of senescence arrest and in cell physiology.

In our related work<sup>40</sup>, some lineages displayed stochastic reversible cell-cycle arrest events, a behavior clearly distinct from type A lineages, which we termed type B. Because we could not reliably propose a characterized molecular mechanism to explain commitment into type B lineages, we decided to take only type A lineages into account for the present study. Type B lineages strongly contributed to experimental
replicative senescence heterogeneity as they represent ~40% of the lineages and showed unpredictable cellular behaviors. However, at the scale of the population of senescent cells, type B lineages are negligible because of their lower fitness. We suggested that they undergo cycles of telomere damage and repair, which could stem from stochastic damage, external stress or mitochondrial dysfunction, globally increasing genome instability. In future works, it would be highly informative to model possible mechanisms of telomere dynamics in type B lineages to quantitatively study their contribution to senescence variability.

In conclusion, it is tempting to speculate that creating a diversity of senescence dynamics for different lineages may be a bet-hedging or a differentiation strategy. Conversely, the heterogeneity generated through telomere shortening may be crucial for genetic, epigenetic or metabolic changes, such as cancer transformation in mammals, which corresponds to rare events occurring to rare cells. Such considerations may not be relevant to wild-type yeast cells as they naturally express telomerase but because replicative senescence is inherently associated with important variability, wired in the structure of telomeres (asymmetry of telomere replication), evolution may have harnessed this variability in multicellular organisms for functional reasons.

## Materials and Methods

**Determination of the asymptotic law of telomere lengths when telomerase is active.** The probability of action of telomerase depends on two parameters $β$ and $L_0$. For our purpose, we want to fit them from only two experimental steady state parameters of telomere length distribution: the mode and the skewness. From a mathematical point of view ($L_n$) is a (time-homogeneous) Markov chain. As the cells have reached a steady state telomere length distribution before telomerase inactivation$^{40}$, we assume that the telomere length distribution is the stationary distribution of this Markov chain. The probabilities of transitions as well as the equation defining the stationary distribution can be explicitly written. By a convergent iterative process, we are able to compute this stationary state. In particular, we can compute its mode and its skewness. By varying the two parameters $β$ and $L_0$, we determine the parameters giving the predefined mode and skewness.

**Determination of the weight $a$ of $L_2$ and of the threshold $Lmin$.** Assume that the parameters $a$ and $Lmin$ are fixed. We pick 24 cells with telomere lengths according to the stationary distribution of the Markov chain ($L_n$). We simulate the evolution of the 32 telomere lengths for each of the 24 subsequent lineages according to the model, and we record for each one the generation of senescence onset. This gives 24 increasingly ordered integers $g_1^l, ..., g_{24}^l$. The integer number $l$ refers to the number of the simulation. We repeat this simulation $N=1000$ or $N=5000$ times.

The error $e(a, Lmin)$ between these numerically obtained generations of senescence onset and the experimental ones, ordered increasingly $G_1, ..., G_{24}$, is measured by the root mean square deviation:

$$e(a, Lmin) = \left[ \sum_{l=1}^{N} \sum_{i=1}^{24} (g_l^l - G_i)^2 \right]^{1/2} \quad (4)$$

We repeat these simulations and calculations for different values of $a$ and $Lmin$.

**Terminal restriction fragment (TRF) Southern blot.** Genomic DNA was extracted from exponential growth cultures using a standard phenol:chloroform:isoamyl (25:24:1) alcohol purification procedure and ethanol precipitation. Then 1 µg of genomic DNA was digested with Xho I and the products were ethanol precipitated, resuspended in loading buffer (10 mM Tris pH 8.0, 1 mM EDTA, 5% glycerol, 0.04% xylene cyanol FF), and resolved on a 1.2% agarose gel for 14 h at 60 V. The gel was then soaked in a denaturation bath (0.4 M NaOH, 1 M NaCl) for 30 min and transferred by capillarity action to a charged nylon membrane (Hybond XL, GE Healthcare). The telomere-specific oligonucleotide probe (5'-GGGTGTGGGTGTGTGTGGTGGG-3') was $^{32}$P-labeled at the 5' terminus with ATP ($\gamma$-$^{32}$P) and T4 polynucleotide kinase (New England Biolabs). The membrane was hybridized using the Rapid-hyb buffer protocol (GE Healthcare). In brief, the membrane was prehybridized at 42°C in Rapid-hyb buffer for 1 h, then the radioactive probe (20 pmol) was added, and the incubation was continued for 1 h. The membrane was washed consecutively with 2 × SSC, 0.5% SDS (42°C for 10 min); 2 × SSC, 0.1% SDS (42°C for 20 min); and 0.1 × SSC, 0.1% SDS (20°C for 30 min). The membrane was then imaged with a Typhoon FLA 9500 scanner (GE Healthcare).

## References

1. McAdams, H. H. & Arkin, A. It's a noisy business! Genetic regulation at the nanomolar scale. *Trends Genet.* **15**, 65–69 (1999).

2. Sumner, E. R. & Avery, S. V. Phenotypic heterogeneity: differential stress resistance among individual cells of the yeast Saccharomyces cerevisiae. *Microbiology* **148**, 345–351 (2002).

3. Thattai, M. & van Oudenaarden, A. Stochastic gene expression in fluctuating environments. *Genetics* **167**, 523–530 (2004).

4. Balaban, N. Q., Merrin, J., Chait, R., Kowalik, L. & Leibler, S. Bacterial persistence as a phenotypic switch. *Science* **305**, 1622–1625 (2004).

5. Nachman, I., Regev, A. & Ramanathan, S. Dissecting timing variability in yeast meiosis. *Cell* **131**, 544–556 (2007).

6. Losick, R. & Desplan, C. Stochasticity and cell fate. *Science* **320**, 65–68 (2008).

7. Süel, G. M., Kulkarni, R. P., Dworkin, J., Garcia-Ojalvo, J. & Elowitz, M. B. Tunability and noise dependence in differentiation dynamics. *Science* **315**, 1716–1719 (2007).
8. Maamar, H., Raj, A. & Dubnau, D. Noise in gene expression determines cell fate in Bacillus subtilis. *Science* **317**, 526–529 (2007).

9. Suda, T., Suda, J. & Ogawa, M. Single-cell origin of mouse hemopoietic colonies expressing multiple lineages in variable combinations. *Proc. Natl. Acad. Sci. USA* **80**, 6689–6693 (1983).

10. Morris, S. A. et al. Origin and formation of the first two distinct cell types of the inner cell mass in the mouse embryo. *Proc. Natl. Acad. Sci. USA* **107**, 6364–6369 (2010).

11. Veening, J. W. et al. Bet-hedging and epigenetic inheritance in bacterial cell development. *Proc. Natl. Acad. Sci. USA* **105**, 4393–4398 (2008).

12. Fujita, M. & Losick, R. Evidence that entry into sporulation in Bacillus subtilis is governed by a gradual increase in the level and activity of the master regulator Spo0A. *Genes Dev.* **19**, 2236–2244 (2005).

13. Spencer, S. L., Gaudet, S., Albeck, J. G., Burke, J. M. & Sorger, P. K. Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. *Nature* **459**, 428–432 (2009).

14. Ansel, J. et al. Cell-to-cell stochastic variation in gene expression is a complex genetic trait. *PLoS Genet* **4**, e1000049 (2008).

15. Raser, J. M. & O'Shea, E. K. Control of stochasticity in eukaryotic gene expression. *Science* **304**, 1811–1814 (2004).

16. Cai, L., Dalai, C. K. & Elowitz, M. B. Frequency-modulated nuclear localization bursts coordinate gene regulation. *Nature* **455**, 485–490 (2008).

17. Smith, J. R. & Whitney, R. G. Intraclonal variation in proliferative potential of human diploid fibroblasts: stochastic mechanism for cellular aging. *Science* **207**, 82–84 (1980).

18. Lundblad, V. & Szostak, J. W. A mutant with a defect in telomere elongation leads to senescence in yeast. *Cell* **57**, 633–643 (1989).

19. Ritchie, K. B., Mallory, J. C. & Petes, T. D. Interactions of TLC1 (which encodes the RNA subunit of telomerase), TEL1, and MEC1 in regulating telomere length in the yeast Saccharomyces cerevisiae. *Molecular & Cellular Biology* **19**, 6065–6075 (1999).

20. Xu, Z., Dao Duc, K., Holcman, D. & Teixeira, M. T. The length of the shortest telomere as the major determinant of the onset of replicative senescence. *Genetics* **194**, 847–857 (2013).

21. Hayflick, L. The Limited *In Vitro* Lifetime Of Human Diploid Cell Strains. *Exp Cell Res* **37**, 614–636 (1965).

22. Proctor, C. J. & Kirkwood, T. B. Modelling telomere shortening and the role of oxidative stress. *Mech. Ageing Dev.* **123**, 351–363 (2002).

23. von Zglinicki, T. Oxidative stress shortens telomeres. *Trends Biochem. Sci.* **27**, 339–344 (2002).

24. Sozou, P. D. & Kirkwood, T. B. A stochastic model of cell replicative senescence based on telomere shortening, oxidative stress, and somatic mutations in nuclear and mitochondrial DNA. *J. Theor. Biol.* **213**, 573–586 (2001).

25. Passos, J. F. et al. Mitochondrial dysfunction accounts for the stochastic heterogeneity in telomere-dependent senescence. *PLoS Biol* **5**, e110 (2007).

26. Sahin, E. et al. Telomere dysfunction induces metabolic and mitochondrial compromise. *Nature* **470**, 359–365 (2011).

27. Epel, E. S. et al. Accelerated telomere shortening in response to life stress. *Proc Natl Acad Sci USA* **101**, 17312–17315 (2004).

28. Blackburn, E. H. Telomere states and cell fates. *Nature* **408**, 53–56 (2000).

29. Proctor, C. J. & Kirkwood, T. B. Modelling cellular senescence as a result of telomere state. *Aging Cell* **2**, 151–157 (2003).

30. Sahin, E. & DePinho, R. A. Axis of ageing: telomeres, p53 and mitochondria. *Nature reviews. Molecular cell biology* **13**, 397–404 (2012).

31. Bahar, R. et al. Increased cell-to-cell variation in gene expression in ageing mouse heart. *Nature* **441**, 1011–1014 (2006).

32. Shampay, J. & Blackburn, E. H. Generation of telomere-length heterogeneity in Saccharomyces cerevisiae. *Proc Natl Acad Sci USA* **85**, 534–538 (1988).

33. Hemann, M. T., Strong, M. A., Hao, L. Y. & Greider, C. W. The shortest telomere, not average telomere length, is critical for cell viability and chromosome stability. *Cell* **107**, 67–77. (2001).

34. Abdallah, P. et al. A two-step model for senescence triggered by a single critically short telomere. *Nat Cell Biol* **11**, 988–993 (2009).

35. Wellinger, R. J., Ethier, K., Labrecque, P. & Zakian, V. A. Evidence for a new step in telomere maintenance. *Cell* **85**, 423–433 (1996).

36. Kelleher, C., Teixeira, M. T., Forstemann, K. & Lingner, J. Telomerase: biochemical considerations for enzyme and substrate. *Trends Biochem Sci* **27**, 572–579 (2002).

37. Olovnikov, A. M. A theory of marginotomy. The incomplete copying of template margin in enzymic synthesis of polynucleotides and biological significance of the phenomenon. *J. Theor. Biol.* **41**, 181–190 (1973).

38. Watson, J. D. Origin of concatemeric T7 DNA. *Nat New Biol* **239**, 197–201 (1972).

39. Lingner, J., Cooper, J. P. & Cech, T. R. Telomerase and DNA end replication: no longer a lagging strand problem? *Science* **269**, 1533–1534 (1995).

40. Xu, Z. et al. Two routes to senescence revealed by real-time analysis of telomerase-negative single lineages. *Nat Commun.* **6**, 7680 (2015).

41. Soudet, J., Jolivet, P. & Teixeira, M. T. Elucidation of the DNA end-replication problem in Saccharomyces cerevisiae. *Mol. Cell* **53**, 954–964 (2014).

42. Teixeira, M. T., Arneric, M., Sperisen, P. & Lingner, J. Telomere Length Homeostasis Is Achieved via a Switch between Telomerase-Extendible and -Nonextendible States. *Cell* **117**, 323–335 (2004).

43. Levy, M. Z., Allsopp, R. C., Futcher, A. B., Greider, C. W. & Harley, C. B. Telomere end-replication problem and cell aging. *Journal of Molecular Biology* **225**, 951–960 (1992).

44. Tan, Z. Intramitotic and intraclonal variation in proliferative potential of human diploid cells: explained by telomere shortening. *J. Theor. Biol.* **198**, 259–268 (1999).

45. op den Buijs, J., van den Bosch, P. P., Musters, M. W. & van Riel, N. A. Mathematical modeling confirms the length-dependency of telomere shortening. *Mech. Ageing Dev.* **125**, 437–444 (2004).

46. Arkus, N. A mathematical model of cellular apoptosis and senescence through the dynamics of telomere loss. *J. Theor. Biol.* **235**, 13–32 (2005).

47. Arino, O., Kimmel, M. & Webb, G. F. Mathematical modeling of the loss of telomere sequences. *J. Theor. Biol.* **177**, 45–57 (1995).

48. Olofsson, P. & Kimmel, M. Stochastic models of telomere shortening. *Math. Biosci.* **158**, 75–92 (1999).

49. Rubelj, I. & Vondracek, Z. Stochastic mechanism of cellular aging—abrupt telomere shortening as a model for stochastic nature of cellular aging. *J. Theor. Biol.* **197**, 425–438 (1999).

50. Dyson, J., Sánchez, E., Villela-Bressan, R. & Webb, G. F. Stabilization of telomeres in nonlinear models of proliferating cell lines. *J. Theor. Biol.* **244**, 400–408 (2007).

51. Portugal, R. D., Land, M. G. & Svaiter, B. F. A computational model for telomere-dependent cell-replicative aging. *BioSyst.* **91**, 262–267 (2008).

52. Olofsson, P. & Bertuch, A. A. Modeling growth and telomere dynamics in Saccharomyces cerevisiae. *J. Theor. Biol.* **263**, 353–359 (2010).
**Acknowledgements**

We thank Erin Henninger for her critical reading of the manuscript. The research of Z.X. and M.T.T. was supported by the European Research Council (ERC-2010-StG 260906—D-END), the Mairie de Paris (Programme Emergences), and the “Initiative d’Excellence” program from the French State (Grant “DYNAMO,” ANR-11-LABX-0011-01). The research of M.D. and T.B. was supported by the European Research Council (ERC Starting Grant SKIPPERAD 306321).

**Author Contributions**

T.B. performed the simulations. M.D. and T.B. designed the mathematical model with inputs from Z.X. and M.T.T. Z.X. performed the experiments and made the figures. M.D. and M.T.T. supervised the work. All authors contributed to the writing of this paper.

**Additional Information**

**Supplementary information** accompanies this paper at http://www.nature.com/srep

**Competing financial interests:** The authors declare no competing financial interests.

**How to cite this article:** Bourgeron, T. et al. The asymmetry of telomere replication contributes to replicative senescence heterogeneity. *Sci. Rep.* 5, 15326; doi: 10.1038/srep15326 (2015).

<figure><img src="image_15.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/