Effects of initial telomere length distribution on senescence onset and heterogeneity.

Replicative senescence, induced by telomere shortening, exhibits considerable asynchrony and heterogeneity, the origins of which remain unclear. Here, we formally study how telomere shortening mechanisms impact on senescence kinetics and define two regimes of senescence, depending on the initial telomere length variance. We provide analytical solutions to the model, highlighting a non-linear relationship between senescence onset and initial telomere length distribution. This study reveals the complexity of the collective behavior of telomeres as they shorten, leading to senescence heterogeneity. Copyright Â

Introduction

Telomeres, the ends of eukaryote chromosomes, are poised in a dynamic equilibrium controlled by two processes: limited telomere shortening at each cell division and elongation by telomerase, a dedicated holoenzyme able to generate de novo telomere sequence. When telomerase is not expressed, as in human somatic cells, or is experimentally mutated in model organisms such as Saccharomyces cerevisiae (Lundblad and Szostak, 1989), telomeres only shorten and after many divisions the cell enters replicative senescence, a permanent cell cycle arrest induced by short telomeres that elicit a DNA damage response. Replicative senescence is implicated in organismal ageing and is a potent barrier to cancer emergence, but its remarkable asynchrony and heterogeneity remain a challenge for investigating the exact relationship between initial telomere length distribution and senescence onset.

Telomere shortening is the unavoidable consequence of the end-replication problem (Olovnikov, 1973, Watson, 1972, Soudet et al., 2014). In most examined species, telomeres end with a 5′ to 3′ singled-stranded DNA overhang (Fig. 1) (Hemann and Greider, 1999, Henderson and Blackburn, 1989, Klobutcher et al., 1981, Makarov et al., 1997, McElligott and Wellinger, 1997, Raices et al., 2008, Riha et al., 2000, Wellinger et al., 1993). When the replication fork reaches the end of the chromosome, the processing of the last Okazaki fragment leaves a gap at the lagging strand, which recreates the single-stranded overhang of the parental telomere (Fig. 1). On the leading strand, after replication, complex maturation steps involving resection and fill-in also regenerate the overhang structure (Larrivée et al., 2004, Faure et al., 2010, Chai et al., 2006, Wu et al., 2012, Soudet et al., 2014). Regardless of these maturation steps, the leading strand template for replication is shorter than the lagging strand one, thus generating after replication two new telomeres of different lengths, one unchanged compared to the parental telomere and the other shorter by exactly the length of the overhang, as illustrated in Fig. 1. Previous mathematical models of telomere shortening also based on the end-replication problem (Levy et al., 1992; Arino et al., 1995; Olofsson and Kimmel, 1999; Arkus, 2005) did not consider the maturation of the leading strand telomere that generates a 3′-end overhang identical to the one on the lagging strand. This maturation step is widely conserved throughout species with the notable exception of angiosperm plants that display a blunt end at the leading telomere (Riha et al., 2000). We also note that other mathematical models examined higher level structures such as t-loops (Griffith et al., 1999, Rodriguez-Brenes and Peskin, 2010), or additional telomere states or breaking mechanisms (Kowald, 1997, Rubelj and Vondracek, 1999, Proctor and Kirkwood, 2002, Proctor and Kirkwood, 2003), such as damage due to oxidative stress (von Zglinicki, 2002). In S. cerevisiae, however, oxidative stress does not significantly alter telomere length (Romano et al., 2013) and the end-replication problem is the main mechanism of telomere shortening.

Fig. 1

Scheme of a chromosome bearing two telomeres and undergoing replication. Telomeres end with a 3′ overhang of length a (measured to be 5–10 nucleotides in yeast (Soudet et al., 2014), chosen here as a=1 or 7 for theoretical or numerical purposes, respectively). After DNA replication, each telomere generates, through either the leading or the lagging strand replication machineries, two new telomeres of different lengths. The coupling effect between the two ends of the same chromosome imposes that only one of the two is shortened while the other retains the parental length.

Consistently, on average, telomeres shorten at a constant rate of exactly half of the overhang length per division. However, while studying the average is informative of the global regulation and homeostasis of telomere length, it misses important contributions of the asymmetry of telomere replication mechanism to the overall telomere length distribution and to the heterogeneity of the onset of senescence. Taking this asymmetry into account, the shortening of a telomere in a cell lineage, defined as a random succession of mitotically related cells (Xu et al., 2015), is probabilistic and follows a Bernoulli process. Additionally, if the two ends of a given chromosome are considered together, the 3′-end at one telomere belongs to the same DNA strand as the 5′-end on the other telomere of the same chromosome, implying that the asymmetry at one telomere is inverted compared to the other (Fig. 1). We define this relationship between the two ends of the same chromosome as a coupling mechanism, which adds another layer of constraint and will also be modeled here.

In this article, we study the consequences of the asymmetry and the coupling on the distribution and the dynamics of telomere length in two distinct phases: at steady state in the presence of telomerase and in a strictly shortening phase without telomerase. We show that the robustness of telomerase recruitment impacts on the variance of the steady-state distribution of telomere length. In turn, this variance defines different regimes of senescence. In a regime of low initial variance, senescence onset cannot be linearly inferred from the average telomere length or even the length of the shortest telomere and we provide an asymptotic expansion to account for this phenomenon. In contrast, a high variance implies a linear correlation between the initial shortest telomere and senescence onset. We provide analytical solutions to the different models we describe and suggest applications for the inference of the initial telomere length distribution from experimental measurements of senescence onset.

Telomeres evolving with telomerase

We first describe the most general model, corresponding to a lineage of haploid yeast cells dividing in the presence of active telomerase. The two telomeres of a given chromosome at generation n, called , are coupled as defined above and the 32 telomeres of the cell shorten according to a Bernoulli random variable of parameter 1/2: if , then L1 is shortened by a nucleotides, whereas L2 is preserved, and conversely if . Telomerase adds new telomere sequences preferentially to shorter telomeres (Teresa Teixeira et al., 2004, Britt-Compton et al., 2009), behavior that we capture by introducing , , Bernoulli random variables of parameter , according to Xu et al. (2013), where f has the shape shown in Fig. 2 (a) and L is the length of the telomere at the extremity i before replication. The shape of f is such that below a length threshold L, the Bernoulli random variable equals 1, that is telomerase is always active. For a telomere longer than L, the probability of decreases to zero, meaning that the longer the telomere, the less likely it is to be elongated by telomerase. Since the number of nucleotides added by telomerase is independent of the length of the telomere (Teresa Teixeira et al., 2004, Xu et al., 2013), we introduce and two independent geometric random variables of parameter p, independent of all the other quantities (including ), which correspond to the number of nucleotides added by telomerase. As a result, for any given chromosome the -valued process follows:where is x if and 0 otherwise. Telomere length is nearly always positive because, in the presence of telomerase, a telomere shorter than is elongated with probability 1, while in the absence of telomerase, senescence is triggered before the short telomere reaches 0 (Bourgeron et al., 2015). Using the same Bernoulli random variable for L1 and L2 mathematically defines the coupling between the two telomeres.

Fig. 2

Steady-state telomere length distribution in the presence of telomerase. (a) and (b) Probability of recruitment and action of telomerase, as modeled from Teixeira et al. (2004) with a length threshold L or simplified with a sharp switch occurring at i. (c) Simulation of the steady-state distribution of telomere length using either (a) (black) or (b) (grey) to describe telomerase recruitment. i was set so as to reach the same mean in the steady-state distribution (Appendix B).

To characterize the steady state of telomere length distribution, we focus on one telomere—because as an approximation, telomeres of different chromosomes are assumed to be independent (Shampay and Blackburn, 1988)—, and consider the projection of the first coordinate of a chromosome in order to compute its equilibrium. We will analyze the coupling effect in more depth in the second regime without telomerase. Our model thus becomes:where L is the length of a given telomere, the length of one of the two daughter telomeres, a geometric random variable of parameter . An averaged version of this model has been studied in Xu et al. (2013); Dao Duc and Holcman (2013) and used in Bourgeron et al. (2015), where instead of being stochastic, telomere shortening was chosen to be deterministic with a constant value of . To make our computations fully explicit without betraying the principles of the biological mechanism, instead of f, we consider a sharp threshold at a value i (Fig. 2b). Our model becomes:

Independently of the value of i, the Markov chain defined by (3) has a unique equilibrium distribution and the generating function of is characterized by the equality:where , and .

Appendix A gives a proof of this result and explains how to get a fully explicit expression for the distribution from Eq. (4). This calculation reveals how the parameters a and p affect the steady-state distribution.

The choice of the value of i according to biological experiments is explained in Appendix B. We find that the variance of the steady-state telomere length distribution obtained using the simplified model (3) is significantly smaller than the one with the complete model (2) (Fig. 2c in black and grey, respectively; 37 bp as compared to 101 bp), demonstrating that the residual recruitment of telomerase to rather long telomeres strongly contributes to the spread of the steady-state distribution of telomere length. In turn, the variance of this distribution is critical for determining the onset of senescence and its heterogeneity, as we show below. Thus, the mode of recruitment and activation of telomerase, dependent on the biochemical properties of the holoenzyme and on its interactions with telomeric proteins (Wellinger and Zakian, 2012), controls key features of senescence once telomerase is removed.

Telomeres evolving without telomerase

We then analyze the consequences of the steady-state distribution on the onset of senescence, meaning the number of generations undergone by a given cell lineage until it enters senescence. We simply call it time of senescence, denoted by T. One practical goal of this section is to derive the parameters of the initial distribution from the time of senescence, which is useful for experimentalists. In senescing cells, telomerase is inactive and when the shortest telomere reaches a threshold S, the cell enters replicative senescence and stops dividing (Abdallah et al., 2009, Armanios et al., 2009, Hemann et al., 2001, Lundblad and Szostak, 1989, Zou et al., 2004). A haploid yeast cell has 16 chromosomes and thus 32 telomeres. Mathematically, we consider the vector of these 32 telomere lengths at generation n. Because each chromosome behaves independently (Shampay and Blackburn, 1988), we can start by studying one chromosome and the behavior of the 16 will easily follow. More precisely, the vector can be seen as a family of 16 independent identically distributed couples each representing the two telomeres of a chromosome, with . The time of senescence is mathematically expressed as:

Normalizations. As telomeres can only shorten, we consider the shortening length to be a=1. For numerical estimations of the time of senescence, we will divide our results by a=7 to obtain biologically relevant values. Moreover, we can choose S=0 by simply translating the initial state by S. These assumptions are made in all following calculations unless stated otherwise.

Distribution of the Time of Senescence. Under these normalization conditions, we find that the distribution of the random variable T is fully explicit:

In particular, its expectation can be written as a function of π as follows:

See Appendix C for the proof. Because of the difficulty to invert this formula, we choose to study separately the influence of the mean and the variance of the initial state on the time of senescence. Thus, we first consider a deterministic and constant initial state :

We define as the first time one of two coupled telomeres reaches zero both starting from x0, and as the time of senescence of the whole cell when the initial state is constant and equals x0.

Almost surely, , this implies that for , and for . For , the law of is given by:

The expected time of senescence is then:

We then perform an asymptotic expansion of for large values of x0, which is numerically justified (Appendix B and Fig. 2c). At the first order the mean behavior prevails:

Concerning the second order, we obtain the following convergence in distribution:See Appendix D for proofs of these results.

The asymptotic development of the time of senescence for one chromosome allows us to derive an approximation of the expected time of senescence (7) by replacing the law of by its asymptotic (9):where erf is the error function defined as:

We find that the expansion (10) is hardly distinguishable from the theoretical process (6), as shown in simulations (compare grey and dashed black lines in Fig. 3 (a)) and can thus be directly used to estimate the mean of the initial state in experimental studies.

Fig. 3

Distinct effects of the mean and variance of the initial distribution on theoretical expressions and numerical simulations of the time of senescence. (a) Starting from a constant distribution , the asymptotic expansion in Eq. (10) is computed and compared to numerical simulations (1000 independent simulations). (b) Starting from a uniform distribution of variance σ and mean , the time of senescence is computed using Eq. (11), which takes only the mean behavior of the initial shortest telomere into account, and compared to numerical simulations (1000 independent simulations).

Influence of the initial variance on the time of senescence. Now, to study only the influence of the initial variance, we consider that each initial telomere is uniformly distributed in the interval and simulate the expected time of senescence as a function of σ (Fig. 3b). When σ has large values, there is a higher probability that the initial shortest telomere of () is far from the mean and, thus, that it remains the shortest one until senescence. We therefore expect, for large enough value of σ (Fig. 3b), that the time of senescence is asymptotically equivalent to the time when the initial shortest telomere reaches zero. As, on average, the number of steps for a simple random walk starting from to reach zero is , we expect the following result:We indeed find this asymptotic behavior by simulations (Fig. 3b), highlighting two regimes that depend on the initial variance. If σ has a small value, which is close to the deterministic initial state studied above (Eq. (10)), the time of senescence is much smaller than expected by just considering the shortest telomere because of the coupling effect. If σ has a large value, the time of senescence is mainly determined by the shortening of the initial shortest telomere.

We next ask which of these two regimes can be observed in simulated times of senescence from the telomere distributions described in Fig. 2, which are biologically more relevant than the previous distribution models. To do so, we simulate 1000 individual lineages of senescing cells and record their time of senescence, starting by randomly drawing their 32 telomeres from the biologically relevant distribution (large variance, in grey in Fig. 2c), from the simplified distribution (intermediate variance, in black in Fig. 2c) or from a constant distribution (no variance) (Fig. 4). We then compare these simulated times of senescence with those predicted either from the mean behavior of the shortest telomere (Eq. (11), dashed black lines in Fig. 4) or the asymptotic expansion on the mean of the initial distribution (Eq. (10), black lines in Fig. 4). The biologically relevant distribution gives simulated times of senescence that are fully predicted by computing the mean behavior of the average initial shortest telomere (Fig. 4a, compare grey and dashed black lines which are superimposed). In contrast, the constant distribution leads to a senescence onset dictated by the asymptotic expansion (Fig. 4c, compare grey and black lines), consistent with the results in Fig. 3a. The simplified distribution produces an intermediate result where the mean behavior of the shortest telomere and the asymptotic expansion lead to similar predictions, which lay close to the simulated times of senescence (Fig. 4b, compare grey, black and dashed black lines). These results show that defining two senescence regimes depending on the initial variance of telomere length distribution is critical for understanding the relevant dynamics of telomere shortening leading to senescence.

Fig. 4

Comparison between simulated times of senescence (grey dots) and predictions from Eqs. (11) (dashed black lines) and (10) (black lines). (a) 32 telomere lengths are randomly drawn from a biologically relevant distribution with a high variance (Fig. 2c, grey distribution) and the time of senescence is simulated to give one data point (grey dot). This process is repeated 1000 times and compared to the two predictions. The grey line represents the average simulated time of senescence. (b) and (c): as in (a), but starting with an intermediate level of variance for the initial telomere length distribution or no variance at all, respectively.

Conclusion

In summary, in this article, we isolated all the sources of fluctuations of the time of senescence that are dependent on telomere length. To do so, we modeled several molecular mechanisms that contribute at various levels to telomere length distribution and dynamics in S. cerevisiae, where they are the most exhaustively and quantitatively described. Among these mechanisms, we found that the asymmetry of telomere replication and the coupling between the two telomeres belonging to the same chromosome significantly contribute to senescence heterogeneity and we formally established their links. We also showed that the mode and robustness of telomerase recruitment control the variance of the steady-state telomere length distribution, which in turn defines two senescence regimes. With a low initial variance, the time of senescence is non-linearly related to the initial mean telomere length. In contrast, a high initial variance leads to a major role of the initial shortest telomere in controlling senescence. Because natural telomere length distributions can vary considerabely, even within a species, we suggest that depending on the initial variance, the two regimes we describe may operate at the same time during senescence. As the core mechanisms modeled here are conserved in most eukaryotes, we expect that our conclusions should also, in principle, apply to telomere-dependent senescence in human cells, although additional factors and mechanisms also contribute to senescence heterogeneity (Griffith et al., 1999, Rodriguez-Brenes and Peskin, 2010, Proctor and Kirkwood, 2002). This work uncovers a new layer of complexity in the relationship between senescence onset and telomere shortening explained by the asymmetry and coupling mechanisms, and proposes methods for assessing the time of senescence or conversely inferring parameters of the initial telomere length distribution.

k	relation
⟦0,a−1⟧	ϕ
a	π_a=21−ppπ_o−∑a−1k=1π_k
⟦a+1,i_s⟧	π_k=−21−ppπ_k−a−1+2−ppπ_k−a
i_s+1	π_i_s+1=−2qπ_i_s−a+(1+q)π_i_s+1−a
⟦i_s+2,i_s+a⟧	π_k=−2qπ_k−a−1+(1+q)π_k−a+qπ_k−1
i_s+a+1	π_i_s+a+1=−2qπ_i_s+π_i_s+1+qπ_i_s+a
>i_s+a+1	π_k=−qπ_k−a−1+π_k−a+qπ_k−1

Table E1

Parameters used in this study.

a	Length of the 3′-end overhang.
S	Threshold length of the shortest telomere inducing senescence.
Ls	Length threshold for the function f, below which telomerase is recruited to the telomere with probability 1.
is	Length threshold for the simplified model of telomerase recruitment.
p	Parameter of the geometric random variable G_n.
β	Parameter of the function f, fitted on experimental data (Xu et al., 2013).