The clock-like accumulation of germline and somatic mutations can arise from the interplay of DNA damage and repair

The rates at which mutations accumulate across human cell types vary. To identify causes of this variation, mutations are often decomposed into a combination of the single-base substitution (SBS) “signatures” observed in germline, soma, and tumors, with the idea that each signature corresponds to one or a small number of underlying mutagenic processes. Two such signatures turn out to be ubiquitous across cell types: SBS signature 1, which consists primarily of transitions at methylated CpG sites thought to be caused by spontaneous deamination, and the more diffuse SBS signature 5, which is of unknown etiology. In cancers, the number of mutations attributed to these 2 signatures accumulates linearly with age of diagnosis, and thus the signatures have been termed “clock-like.” To better understand this clock-like behavior, we develop a mathematical model that includes DNA replication errors, unrepaired damage, and damage repaired incorrectly. We show that mutational signatures can exhibit clock-like behavior because cell divisions occur at a constant rate and/or because damage rates remain constant over time, and that these distinct sources can be teased apart by comparing cell lineages that divide at different rates. With this goal in mind, we analyze the rate of accumulation of mutations in multiple cell types, including soma as well as male and female germline. We find no detectable increase in SBS signature 1 mutations in neurons and only a very weak increase in mutations assigned to the female germline, but a significant increase with time in rapidly dividing cells, suggesting that SBS signature 1 is driven by rounds of DNA replication occurring at a relatively fixed rate. In contrast, SBS signature 5 increases with time in all cell types, including postmitotic ones, indicating that it accumulates independently of cell divisions; this observation points to errors in DNA repair as the key underlying mechanism. Thus, the two “clock-like” signatures observed across cell types likely have distinct origins, one set by rates of cell division, the other by damage rates.


Introduction
Mutations are the net result of multiple processes, including endogenous and exogenous DNA damage, errors in DNA replication, and DNA repair.The rate at which they accumulate varies substantially across the human body: mutation rates in germline lineages are typically very low, averaging fewer than one mutation per haploid genome per year, whereas hundreds of mutations accrue yearly in tissues exposed to extensive exogenous mutagens [1][2][3].The origins of these pronounced differences remain poorly understood: while the biochemical mechanisms that underlie mutations have been well characterized (e.g., [4,5]), their relative importance in a given cell type or tissue remains largely unknown.
Over the past decade, a new approach to this question has become feasible, with the analysis of whole genome mutational spectra in cancer genomes and the identification of repeatable patterns of mutations termed mutational signatures [6][7][8].In these analyses, single base substitutions (SBSs) are classified into 96 types based on the identity of the substitution and the flanking nucleotides, and each tumor sample is characterized by the distribution over the 96 types.Variation across samples is then modeled as a linear combination of relatively few signatures, each characterized by a different profile of the 96 SBS types.The proportion of mutations attributed to a given signature is often described as its "activity" or "exposure" (e.g., [8,9]), though in what follows, we use the term loading.The idea behind this decomposition is that the signatures reflect one or more mutational processes occurring in different types of samples.In practice, signature profiles and signature loadings across samples are inferred jointly using nonnegative matrix factorization [6].Nearly 100 signatures have been identified to date [8], which are cataloged and updated in the COSMIC database [10].Only a subset of COSMIC signatures have a known or partially known etiology, while the sources of most remain elusive.Part of the difficulty is that all signatures reflect the interplay between the initial DNA damage or replication error and DNA repair, and disentangling these from their net effects is not straightforward.
In studies conducted to date, two such COSMIC SBS signatures, signature 1 (SBS1) and signature 5 (SBS5), are ubiquitous.They are detected in tumor samples as well as in healthy somatic cells in humans and other mammals [11], in many cell types contributing the majority of mutations [2].They are also the predominant contributors to the spectrum of germline mutations [12].The number of mutations attributed to SBS1 correlates with the age of cancer diagnosis in many cancer types, as do mutations assigned to SBS5, leading the 2 signatures to be described as "clock-like" [13].Because they share this behavior, mutations attributed to these 2 signatures are sometimes combined and analyzed together [14][15][16].
SBS1 is dominated by cytosine transitions in the CpG context.The rate of SBS1 accumulation varies among cell types and increases with the rate of cellular turnover [3,17].The loading of SBS1 mutations has also been found to increase in metastatic tumors relative to their primary tumor counterparts, plausibly due to accelerated cell division rates in metastasis [18].
In contrast, SBS5 has a more diffuse distribution among the 96 SBS types.Its ubiquity across cell types, cancerous and healthy, has led to the suggestion that the signature reflects "background" endogenous cellular damage to DNA [19,20].However, the number of SBS5 mutations is also affected by the presence of at least one exogenous mutagen, tobacco smoke [21,22].
Here, we consider conditions under which mutations can accumulate in a clock-like fashion and how the rate of accumulation depends on underlying processes of damage and repair as well as DNA replication.Analyzing whole-genome sequencing data sets of somatic and de novo germline mutations in cell types with varying rates of cell division, ranging from postmitotic neurons to rapidly dividing intestinal epithelium, we identify mechanisms that likely underlie the clock-like mutation accumulation of SBS1 and SBS5 and highlight key differences between them.

General model of mutation accumulation
We study the interplay between DNA replication, damage and repair, by extending the model introduced in [23].We first consider the mutational processes that occur between cell divisions (Fig 1A).A given source of damage causes lesions at a rate u per base pair, which are detected and repaired at a rate r.Repair leads to mismatches with probability �, due to the misincorporation of a nucleotide.We assume that once repair is complete, there are no mechanisms that differentiate the newly synthesized strand from the strand used as a template by the DNA repair polymerase.We denote the rate of mismatch resolution by q and the probability of correct resolution by p.The outcome of mismatch resolution can vary depending on the type of the mismatch and the local sequence context.We discuss the special case of T:G mismatches in CpG contexts in more detail below.
Based on these assumptions, we derive the expected number of lesions, mismatches, and mutations and their variances as a function of the time since the last round of DNA replication t in a genome of length l.The full analysis of the model is presented in Methods; here, we describe the main findings.Over short time periods, t<r −1 , during which repair has had little chance to occur, the number of lesions grows at rate ul and the number of mismatches at rate ul�.Over longer times, t�r −1 and �q −1 , the expected number of lesions and mismatches reach a steady state, at which they are approximately equal to ul/r and ul�/q, respectively.In turn, the number of mutations due to incorrect repair increases at rate ul�(1−p).
From independent lines of evidence, we know that the vast majority of DNA damage in healthy cells does not lead to mutations: the numbers of mutations in healthy cells are substantially lower than estimated damage rates would suggest [24], and mutation rates in individuals with DNA repair deficiencies are orders of magnitude higher [25].These observations imply that in healthy cells, the repair rate is on the order of the total damage rate, i.e., that r*ul.Therefore, the steady state between damage and repair likely is established long before the next cell division.For simplicity, we further assume the rate of mismatch resolution is of the order of the repair rate, q*r.We note that the dynamics of cancer cells may be different from those of the healthy cells on which we focus here; in particular, the repair machinery may be debilitated or overwhelmed and lesions may persist over multiple cell divisions, a phenomenon termed "lesion segregation" [26,27].
The next part of the model describes the mutational processes that occur during replication (Fig 1B).For simplicity, we assume that the cell divides immediately after DNA replication and treat the two processes as simultaneous, occurring at a fixed rate ϕ.An unrepaired lesion stalls DNA replication and triggers the recruitment of polymerases that can replicate over the lesion, through translesion synthesis [28].This synthesis leads to the incorporation of the incorrect nucleotide opposite the lesion on the template strand with probability R, which depends on the type of lesion and possibly on replication timing [29].If, with probability 1-R, the correct nucleotide is incorporated then, given our assumption that repair is rapid relative to the rate of cell division (i.e., that r�ϕ), the lesion is likely to be repaired during the next cell cycle.If the translesion polymerase incorporates an incorrect nucleotide, however, the repair process will propagate the error to the complementary strand, generating a mutation.Overall, erroneous translesion synthesis introduces mutations at a rate of ulRϕ/2r.
Mismatches unresolved during the cell cycle cause a mutation in one of the 2 daughter cells.Given our assumption that the rate of mismatch resolution far exceeds the rate of cell division (i.e., q�ϕ), these mutations track cell division and accumulate at a rate of ul�ϕ/2q.In contrast, mutations caused by repair errors are independent of cell divisions and accumulate with absolute time, at a rate of ul� (1−p).
Lastly, mismatches also arise during DNA replication due to the misincorporation of nucleotides by replicative polymerases.We denote the probability of a misincorporation per base pair by w and the probability that such a misincorporation leads to a mutation by P. Importantly, replicating DNA carries transient features that distinguish the newly synthesized strand from its template [30] and help mismatch repair substantially decrease the number of replication errors that become mutations, i.e., P�1.In sum, the number of mutations due to replication errors increases with age at a rate of wlPϕ.
Considering these different processes together, our model predicts that the expected number of mutations m at age x is given by The two terms in this equation correspond to two kinds of clock-like behaviors.The first depends on the cell division rates and includes damage-induced mutations as well as DNA replication errors; assuming, as we do, that cell division rates are fixed, these mutations accrue with age at a constant rate.The second type of mutation is driven by errors of DNA repair in response to damage; assuming damage rates are constant, it too depends on absolute time.However, a distinguishing feature of the two types of clock-like mutations is how they behave as a function of cellular turnover rates; in particular, postmitotic cells should show no increase of the first type of mutations.
Accumulation of CpG transitions.The process of spontaneous cytosine deamination contributes substantially to the mutation rate [31].The standard explanation is that because at methylated cytosines, the deamination rate results in a thymine, one of the canonical bases, the efficiency of repair is low [32].To investigate the dynamics underlying the accumulation of CpG transitions, we first consider the consequences of methylated cytosine deamination during the cell cycle.Given that deamination leads directly to a mismatch, we can employ a simpler model, originally introduced in [23].The dynamics is analogous to the general model and can be recovered in the q!0 limit (see Methods, Eq 12).
Methylated cytosines in double-stranded DNA deaminate spontaneously at a high rate (estimated as u d = 2.3×10 −5 per year [33]).The resulting T:G mismatches can be detected and repaired by base excision repair (BER, [34]) or mismatch excision repair (MMR, [35]) at effective rate r.With probability �, the repair machinery erroneously substitutes the guanine for an adenine, leading to a mutation (Fig 2A).
As in the more general model, at early times, t<r −1 , the number of mismatches grows linearly with the rate u d , eventually reaching a steady state with the expected number of mismatches equal to u d l/r, where l stands for the total number of methylated cytosines.Given the high estimated rate of spontaneous deamination and the relatively low observed rates of SBS1 mutations [13], repair must be efficient and have time to act before the cell division.The implication is that the repair rate is higher than the cell division rate, r�ϕ, and that the cell divides after the number of mismatches has reached a steady state.At DNA replication, unrepaired mismatches lead to mutations in one of the two daughter cells.Therefore, the number of mutations m at age x depends on the cell division rate ϕ.On the other hand, the number of mutations due to defective repair follows absolute time, and such mutations accumulate regardless of cell divisions.
Under an alternative hypothesis that is not mutually exclusive (Fig 2B ), CpG transitions result from the deamination of methylated cytosines in single-stranded DNA immediately prior to replication.Mismatches arising at this point cannot be repaired.The number of mutations that arise from such deamination events is proportional to the expected number of cell divisions ϕx and the time of transient single-strandedness Δt.The rate of deamination of methylated cytosines in single-stranded DNA is estimated to be orders of magnitude higher than in double-stranded DNA (u s = 3.5×10 −3 per year [33]), making this scenario a plausible explanation for the observed SBS1 mutation rates: to account for 1 mutation in 10 divisions, the time of transient single-strandedness would need to be of order Δt*1 min.This order of magnitude can be compared with the typical velocity of replication fork in human cells, of order 1 kb per minute [36].
Some of the CpG transitions could arise from DNA replication errors, if methylated cytosines are a difficult template for DNA polymerases [37,38].On average, replicatory errors contribute wlP transitions per cell division, where w denotes the probability of a misincorporation per base pair, and P the probability that such a misincorporation is left unrepaired by mismatch repair and leads to a mutation.
Taken together, the expected number of CpG transitions m at age x is given by where l denotes the number of methylated cytosines in a CpG context.
Relating model predictions to data.In any given cell, a combination of clock-like and non-clock-like mutagenic processes contribute to the mutation rate.Thus, although in our model, we focus on the conditions under which a clock-like accumulation of mutations will arise, in practice, the total number of mutations of any signature observed in a cell will likely not be strictly proportional to time (Eqs 1 and 2).Instead, there may be a non-negligible contribution of mutations that do not increase with age (i.e., are "age-independent").These mutations can have several different sources: for example, they may have accrued during the first few cell divisions of development [39], the last stages of cell differentiation [21], or in response to an acute exposure to exogenous mutagens [40].If these mutations occur in a burst, i.e., over a short period of time, then they will contribute a set number of mutations regardless of age, and a regression of the number of mutations against age can lead to a positive intercept.
A nonzero intercept can also be generated by a clock-like process occurring at a varying pace over development.We discuss model expectations in this case with a toy piece-wise linear model, in which mutations accumulate at 2 different rates over 2 time periods (S1 Fig):  [23].Double-stranded cytosines deaminate at a constant rate u d per base pair and the resulting mismatches are repaired at a rate r.With probability �, the mismatch resolution is incorrect and leads to a mutation.We assume that the cell divides immediately after DNA replication and treat the two processes as simultaneous (occurring at a rate ϕ).Unresolved mismatches at cell division lead to a mutation in one of the two daughter cells.We provide the prediction of the model for the number of mutations m at a given age x.Inefficiently repaired mismatches accumulate with the cell divisions, which occur at rate ϕ, and repair errors accumulate with absolute time, independent of cell divisions.The number of methylated cytosines is denoted by l. (B) An alternative source of CpG transitions is the deamination of methylated cytosines in single-stranded DNA (at a rate of u s ) during DNA replication.The model assumes that a deamination immediately preceding the polymerization of the second strand is not repaired and leads to a mutation in one of the daughter cells.The expected number of mutations is proportional to the rate of cell divisions and the time of transient single-strandedness Δt.
https://doi.org/10.1371/journal.pbio.3002678.g002initially, the number of mutations grows at rate μ 0 , and after some time x 0 , the underlying parameters of mutagenesis (Eq 1) change and mutations accumulate at rate μ6 ¼μ 0 .For data points collected after time x 0 , the expected number of mutations is then where a = μ and b = (μ 0 −μ)x 0 .Hence, the intercept is positive if μ 0 >μ; for example, if entering a postmitotic state at x 0 .
Conversely, an increase in the rate of cell division or the rate of damage at x>x 0 (such as occurs in the lung epithelia of regular smokers [21]) will lead to an increased slope in the number of mutations with age.If the mutation rate was significantly lower at earlier stages (say, before the person smoked), i.e., if μ 0 <μ, then a regression of the total number of mutations on age will yield a negative intercept.
In practice, another reason for a positive intercept in data are sequencing errors that arise from technical artifacts and will not depend on the age of the donor.
In reality, more than one of these phenomena is likely operating at once, and it is therefore important, in disentangling the origins of clock-like mutations, to distinguish those that accrue with age from those that contribute to the intercept, and may or may not be clock-like.Concretely, in data sets for which there is a significant positive intercept, the goal is to tease apart signatures that contribute at a constant level in all samples from signatures for which mutations increase in number with age.To this end, we extend the standard signature decomposition method to allow for a mixture of 2 components, a constant and an age-dependent one, and attribute mutational signatures to the two jointly.

Clock-like signatures across cell types
In order to gain insight into the origins of mutations that accumulate with age, we analyze patterns of mutation accumulation across cell types with different characteristics.To this end, we consider data sets that provide single-cell resolution mutation data, collected using a variety of experimental approaches (see Methods), including mutations in neurons and muscle cells [3], liver hepatocytes [41], lung epithelium [21], small bowel epithelium [42], colonic epithelium, and testis seminiferous tubules [2], as well as germline mutations identified from blood samples of pedigrees [43,44].We rely on mutation data from donors without a disease diagnosis and on lung samples from non-smokers (see Methods).
For each cell or tissue type, we attribute mutational signatures by relying on the COSMIC database of signatures inferred from a large collection of cancer samples [8], as also done previously to describe mutational landscapes in noncancerous soma (e.g., in ref. [41]) and germline mutations (e.g., in [12]).Most of these signatures have been linked to specific mechanisms or are associated with exposures to mutagens, and they therefore provide a useful basis for analyzing and comparing mutation accumulation across tissues and cell types.Nonetheless, because the signatures were originally inferred from tumor samples, they may not fully capture the mutational processes acting in normal cell types, particularly in the male and female germ cells.This limitation could lead to a poorer fit, as well as to incorrect assignments of mutations to signatures to which they do not, in fact, belong.Here, we focus on clock-like signatures and choose a method of signature attribution that limits erroneous assignments of mutations to SBS1 and SBS5 and thus avoids overestimating their contributions (see Methods for details).
To compare observations with our model predictions, we develop an approach to focus on the subset of mutations in a given cell type that accumulate with age.This is done in 2 steps.First, we fit a linear model for all mutations jointly, i.e., y = ax+b, where y denotes the number of mutations per genome and x denotes age (see Fig 3A for mutations in neurons).Second, we model the distribution over the 96 substitution types as a mixture distribution of 2 components: the constant component (Fig 3B , yellow) contributes on average the same number of mutations in samples of all ages, whereas the age-dependent component (Fig 3B , blue) contributes an increasing number of mutations with age.We decompose the slope and the intercept into COSMIC signatures jointly, such that the dependence of the number of mutations y S attributed to a given signature s on age takes the form where P a (s) and P b (s) denote the loadings of signature s in age-dependent and constant components, respectively.We estimate the loadings by extending the standard methods of signature attribution (see Methods); see Fig 3C and 3D for the example of loadings estimated for mutations in neurons.We note that applying this method is only possible if the intercept b is large enough (i.e., if the data contains enough age-independent mutations to attribute mutations to signatures in the constant component).In all cell types except for colonic epithelium, we find significant positive intercepts in the regression of the number of mutations on age, consistent with a burst of mutations in early development, for example.In these cases, we decompose mutations into the age-dependent distribution P a and the constant component distribution P b .In the case of colonic epithelium (Fig 3L ), the intercept is significantly negative, possibly because the mutation rate during ontogenesis is lower than in adult life (S1B Fig) , when the exposure to damage and the cell division rate is higher [45].To proceed with our analysis in this case, we assume (rather than infer, as for other cell types) that all signatures increase at constant rates with age.
Our analysis reveals multiple signatures that increase with age (see S2 Fig) , including the 2 that were previously reported [13], SBS1 and SBS5, and 2 additional signatures that are common across cell types considered, SBS12 and SBS16.SBS16 may not be independent from SBS5 [8,10].In turn, SBS12 contributes up to 7% to 10% of mutations in neurons and the female germline, approximately 7% of mutations in lung and liver, approximately 3% in small bowel and colon, approximately 1% in the paternal germline, but is not found at detectable levels in muscle.Both SBS12 and SBS16 are of unknown etiology and dominated by T to C/A to G transitions.
In examining the mutation accumulation across cell types that vary in their division rates, we focus on SBS1 and SBS5, the two ubiquitous clock-like signatures, and use our decomposition method to examine possible sources for their age dependencies.If driven by cell division, we predict that the rate at which they will accumulate should vary substantially with cellular turnover rates.In contrast, if driven by damage rates, the rate should be much less sensitive to turnover rates, but may vary among tissues owing to differences in endogenous and exogenous damage rates.
In this regard, the accumulation of mutations with age observed in neurons is particularly informative, given that neurons are fully postmitotic cells.Despite the lack of cell divisions, mutations accumulate at rates similar to actively dividing lineages [3].Using our decomposition, signatures whose mutation numbers increase with age are distinct from those that do not (Fig 3C and 3D).Notably, the increase with age is predominantly driven by mutations assigned to signature SBS5, with secondary contributions from SBS16 and SBS12.Strikingly, there is no discernible contribution of SBS1, as we discuss in more detail below.In turn, mutations in the constant component are attributed primarily to signatures SBS5 and SBS1, as well as signature SBS89, which is of unknown etiology but has been reported to be active in the first decade of life [45].
Mutations assigned to the clock-like signature SBS5 are found across cell types and increase significantly with age in every one (Fig 3E -3L).Moreover, SBS5 is the prevalent mutation signature in all cell types, except for small bowel and colon, for which more mutations are attributed to SBS1.That SBS5 is the dominant signature in postmitotic cells such as neurons, as well as in maternal mutations, most of which arose in oocytes, indicates that such mutations can arise independently of DNA replication cycles and points to errors in DNA repair, which accumulate with damage rates (Eq 1).Similarly, SBS12 and SBS16 contribute to both neurons and female germline mutations as well as to mutations in rapidly dividing cells (S2 Fig) , suggesting that the age dependencies of these signatures are not driven by DNA replication cycles either.
When not arising from replication errors, our model predicts that the number of mutations will be clock-like only if the damage rate u is constant.If we assume that probabilities � and p are fixed, as seems sensible if they are primarily determined by inherent properties of DNA repair (e.g., the error rate of a polymerase), then the variation in the rate of SBS5 mutation across cell types reflects differences in rates of endogenous and exogenous damage.Consistent with this notion, the rate of SBS5 mutation is highest for epithelia in the colon and lung (Fig 3), which plausibly experience high rates of damage, and lowest for mutations assigned to the maternal genome, potentially reflecting the fact that oocytes are particularly well protected [46,47].This model also helps to explain the observation that increasing the damage rate by exogenous factors, such as long-term exposure to tobacco, significantly increases SBS5 mutation rate in lung cells [21,48].
In that light, it may seem puzzling that in such different cell types, which presumably experience distinct sources of damage, a large fraction of mutations are consistently comprised of SBS5.As an explanation, we propose that SBS5 reflects errors in DNA synthesis during repair, a critical step in many repair pathways (e.g., nucleotide excision repair or homologous recombination) [4].These pathways often involve the synthesis of multiple nucleotides surrounding the lesion, using the intact strand as a template.The errors of the gap-filling polymerase may be displaced from the position of the original lesion, disassociating the mutational signature from the context of the original damage.We therefore hypothesize that this mutational signature reflects the error profile of the polymerase (�) and the asymmetry of mismatch resolution (p).
The second signature to increase with age, SBS1, does so in all cell types considered, except for liver, where hepatocytes are routinely dormant in the cell cycle [49], and neurons, which are postmitotic.More generally, the rate at which mutations assigned to SBS1 increase with age varies widely among cell types and is highest in those characterized by the highest turnover rates (such as intestinal epithelia, where turnover time estimates are of the order of 3 days [50]).Thus, SBS1 appears to be driven by cell division rates.A possible exception is the observation of a slight increase with age in maternal germline mutations, most of which arose in oocytes (see the discussion of germline mutations below).These observations are in agreement with previous observations from cancer studies [6,13,51].
The origin of CpG transitions remains unclear.If they arise because methylated CpGs are a poor template for replication or from spontaneous deamination of single strands during replication, their dependence on cell divisions is expected.Less intuitively perhaps, the same expectation holds if they arise from spontaneous deamination and are efficiently and accurately repaired during the cell cycle (Fig 2).Current data do not allow us to pinpoint when in the cell cycle the damage accrues, however.Two plausible sources are unrepaired mismatches that accrue during the cell cycle and deamination of single-stranded cytosines during DNA replication.As we show, their relative importance will depend on efficiency of mismatch repair as well as the length of time spent single-stranded during replication, parameters that are to our knowledge unknown.
Regardless, we can use the fact that SBS1 does not discernibly increase with age in postmitotic neurons in order to estimate an upper bound on the rate of SBS1 mutations due to repair errors.Given no cell divisions, ϕ = 0, we estimate the upper limit of the error rate of the resolution of T:G mismatches at CpG sites to be � � 5% m x u d lÞ (here, we assume the detection threshold to be �5%, the lowest loading of an attributed signature in neurons).Transient single-strandedness, as could possibly arise during transcription, or double strand break repair [52], could enhance the rate of deamination [33], in which case the error rate of repair would need to be lower.Our estimate is similar to the measurements of the fidelity of polymerase beta [53], employed by BER, one of the pathways that repairs T:G mismatches [34].It is also of the order of the lower bound on the error rate for DNA synthesis without proofreading, � 0 *10 −4 , found by considering the equilibrium kinetics of DNA synthesis, given the energy difference between a mismatch and a correct base pair [54,55].These calculations show that, despite a high rate of spontaneous deamination, repair mechanisms should be accurate enough for incorrect repair to be an insignificant source of mutations.Instead, repair of this type of mutation is likely both very efficient and accurate in all cell types, leading the number of SBS1 mutations to track cell divisions.

Mutation accumulation in the germline
While SBS1 and SBS5 are known to predominate among germline mutations [12,40], we expect the rates of mutation accumulation to differ between the sexes, given the pronounced differences in gametogenesis.In mothers, any mutation that increases with maternal age should have arisen in an oocyte, a postmitotic cell (although a small fraction may also arise in the early development of the child, if children of older mothers have more mutations in the first few cell divisions [56]), whereas in fathers, age-dependent mutations should arise in dividing spermatogonia.In turn, the mutations that do not depend on parental ages in either sex originated either prior to the onset of puberty in the parents or soon after fertilization of the offspring [57].Given that in humans, early development is the same in both sexes until the *6th week of embryonic life [58], we might expect some similarity between the mutation types that contribute to the constant component in the two sexes.
In agreement with these expectations, for maternal mutations, the distributions of mutational signatures differ markedly between the age-dependent and constant components: in the constant component distribution (Fig 4D and 4H), C to T/G to A transitions dominate, with leading contributions of SBS1, SBS6 (associated with defective DNA mismatch repair [8]), and SBS30 (associated with defective base excision repair [59]).The age-dependent distribution is significantly more diffuse (Fig 4C and 4G), with top contributions from SBS5, SBS12, SBS16, and SBS39.SBS39 predominantly features C to G/G to C transversions, a substitution type known to increase sharply with maternal age and associated with double strand break repair [56,60].Given that mutational signatures have been identified from cancer somatic tissues [8], it is conceivable that SBS39 absorbs a significant portion of the C to G/G to C substitutions characteristic of maternal mutations, even if the process that generated them in the germline has a distinct etiology.Qualitative conclusions are similar when adding the father's age as a covariate in the model (see S3 Fig).
Surprisingly, there is a small but significant increase of SBS1 mutations with maternal age (0.019 mutations per year per gamete), which seems at odds with the lack of cell divisions in oocytes.One possibility is that there is an increase in the steady state number of T:G mismatches in aging oocytes, potentially due to reductions in repair efficiency r in older mothers [56].Uncorrected T:G mismatches in oocytes would lead to zygotic mutations in the first division after fertilization and be detected by pedigree sequencing of trios.Alternatively, the slight increase of SBS1 mutations with age in maternal mutations, in contrast to what is seen in neurons, could be explained by a higher error rate of the repair of T:G mismatches in oocytes.In other words, while the rate of erroneous repair of mismatches due to cytosine deamination may be minimal in neurons, it could be higher and detectable in oocytes.In principle, a similar outcome would be observed if the rate of erroneous repair is the same but cytosines deaminate more often in oocytes relative to neurons, but this explanation seems less likely given the lower levels of genome-wide DNA methylation in oocytes until they enter the growth phase [61].
Among paternal mutations, the age-dependent and age-independent distributions are both dominated by signature SBS5.Most of the mutations in the constant component (Fig 4F) are attributed to SBS5, but there is also an enrichment of C to T/G to A transitions (SBS32, SBS1, and SBS30).Unlike in mothers, the increase with age is driven by SBS5 and SBS1 jointly (Fig 4E).To test whether age-dependent paternal germline mutations accumulate primarily in spermatogonia, as expected, we compare these findings with the decomposition of somatic mutations detected in the seminiferous tubules of the testes [2].As reported by [2], the slope of the regression line is very similar for the 2 data sets (Fig 4B).In contrast to their study, however, the intercept for testis is not significantly different from 0, while the intercept for germline mutations is of approximately 10 mutations per haploid genome; the reasons for this difference are unclear to us.Regardless, in our analysis, the age-dependent component of paternal mutations has a very similar distribution of signatures to those inferred for testis (Fig 4E and inset).
Overall, the age-dependent distribution is highly similar in the two sexes (Fig 4G and 4I), except for SBS1, which is significantly more pronounced in the paternal germline.The similarity of age-dependent mutational signatures between the two sexes suggests that it is not only the continuous cell divisions but also a higher damage rate that distinguishes the paternal germline from the maternal one.While the constant components differ between maternal and paternal mutations, SBS1 and SBS30 are found in both, likely reflecting a contribution of gonosomal mutations.

Discussion
As we show by modeling, distinct mutagenic processes can give rise to the clock-like accumulation of mutations observed in the germline and soma, so long as cell divisions and damage occur at a reasonably constant rates.To tease these processes apart, we estimate the rate of accrual of clock-like mutations across dividing and non-dividing cell types.Our analysis reveals that the ubiquitous SBS1 and SBS5 originate predominantly from different sources: whereas SBS1 tracks cell divisions, SBS5 accumulates in postmitotic as well as dividing cells, and appears to track DNA damage levels.
Based on the behavior of SBS5 across cell types, we hypothesize that such mutations arise from errors in DNA synthesis during repair.This hypothesis could be explored further by examining, for example, if SBS5 mutations are enriched at loci with high repair rates.In turn, the dependence of SBS1 on cell divisions suggests that the rate of accumulation of such mutations could serve as a "counter" for cell divisions [13,62], applicable to different cell lineages.For now, however, cell division rates at different stages of human development remain poorly characterized, stymieing such efforts.
In gaining a better understanding of cell division rates, it will be interesting to examine if rates of DNA damage and cell division covary.As an example, reactive oxygen species-an important source of DNA damage-influence cell cycle progression [63].More broadly, cell metabolism and cell cycle progression are tightly related (reviewed in [64]).This relationship could contribute to the observed variations in clock-like mutation rates among different cell types within an organism and potentially to differences across species.For example, cells in smaller, shorter-lived mammals typically exhibit higher mass-specific metabolic rates [65], shorter cell cycles [66], and higher mutation rates for both SBS1 and SBS5 [11].While establishing causality remains a challenge, these findings hint at a potential link between rates of DNA damage and of cell division.
Our analyses further confirm the importance of damage as a source not only of somatic mutations, but also for germline mutations [56,67,68]: over two-thirds of germline mutations are assigned to SBS1 and SBS5, signatures that arise from damage that is either not repaired or repaired incorrectly.The source of such DNA damage is most likely endogenous cellular processes, accounting for the omnipresence of these two signatures across cell types [2] and species [69], as well as their characteristic clock-like behavior under most conditions [48].
While our analyses help to make sense of a number of disparate observations in human germline and soma, they also raise a number of new questions.In particular, it is puzzling that CpG transitions accumulate with absolute time in phylogenetic data from mammals [70,71], when SBS1 mutations accumulate with cell divisions within humans, and there are dramatic differences in germ cell division rates across mammalian species.
It is also unclear how species-specific rates are set for other types of mutations.The dominance of SBS5 in the germline implies that most of the variation in mutation rates across species is likely explained by differences in the rates of DNA damage and repair (e.g., [72,73]).To what extent this balance is directly shaped by selection versus a byproduct of changes in cellular activity remains to be explored, however.

Model of mutation accumulation
We first consider the interplay of damage and repair during the cell cycle.The time since the last DNA replication is denoted by t.We denote the number of intact sites by n 0 , the number of lesions by n 1 , the number of mismatches by n 2 , and the number of mutations by n 3 , where all these sum up to the genome size, i.e., n 0 +n 1 +n 2 +n 3 = l.DNA damage leads to new lesions at rate u.Lesions are repaired at rate r, with error rate �.Incorrectly repaired lesions lead to mismatches, which are resolved at rate q.We assume the mismatch resolution mechanism cannot discern the correct base in the mismatch site and it results in a correct base pair with probability p (see schema in Fig 5).The model parameters and variables are summarized in Table 1.
The expected number of sites in each state is well approximated by the following system of equations: where we have neglected the probability of damage affecting the same site multiple times, as is realistic.The number of intact sites is always much greater than the number of lesions,  mismatches, and mutations combined, i.e., n 0 �l�n 1 +n 2 +n 3 .This dynamics is therefore well approximated by The solution to this system with initial condition n 0 (0) = l (and thus n 1 (0) = n 2 (0) = n 3 (0) = 0) is At steady state, which is achieved at time t�r −1 and q −1 , the expected numbers of lesions and mismatches are constant, and the number of mutations n 3 grows linearly with time.Specifically, where rþq rq is the delay due to lesion and mismatch processing.We assume, as is plausible, that a steady state is established before cell division, which in terms of model parameters means that the rate of cell division obeys ϕ�r and q, and consequently rþq rq is negligible at timescales of ϕ −1 or longer.
Next, we consider the expected number of mutations m over a period x that spans multiple cell divisions occurring at rate ϕ.We additionally account for unrepaired replicatory errors, of which there are wlP per cell division, where w denotes the error rate of a replicative polymerase and P is the probability that a mismatch is unrepaired by mismatch repair.Altogether, the expected number of mutations is given by where we assume that unrepaired lesion leads to a mutation in one of the two dividing cells with probability R.
To quantify variation in the number of mutations, we note that the set of Eq (6) corresponds to a monomolecular reaction system and therefore that the equivalent chemical master equation is solved by a product Poisson distribution [74].Thus, the variances of the random variables N 1 ,N 2 , and N 3 (numbers of lesions, mismatches, and mutations) are equal to their expected values and given by We further assume that cell division dynamics is independent of the dynamics within a cell cycle, and that the number of cell divisions in time x is Poisson distributed with mean ϕx.We can therefore estimate the variance of the number of mutations M at age x using the law of total variance, To obtain the reduced model describing the dynamics underlying CpG transitions, we take the limit q!0.This results in where n 1 now stands for the number of T:G mismatches, and n 2 is the number of mutations.This system is solved by Unrepaired mismatches lead to a mutation in one of the 2 dividing cells.Additionally, we consider that cytosines can deaminate at rate u s during the transient single-strandedness that occurs at DNA replication, for time Δt.Analogously to the general model, assuming the dynamics reaches a steady state before the cell division, the expected number of mutation is then given by

Data analysis
We analyze the patterns of mutation accumulation in different cell types, using a variety of publicly available data sets that each include 7 or more individuals of different ages.The studies used a variety of approaches in order to characterize mutations that accrue in a given cell type, including single molecule sequencing (S), sequencing of small monoclonal samples (C), and sequencing of cell colonies derived from a single cell (D).In addition, we consider de novo mutation calls based on sequenced genomes from blood samples of human trios (P).In this case, mutations are assigned to maternal and paternal genomes using "read-backed phasing" [60,75] as well as by transmission to a third generation [44,60]; such mutations are assumed to have arisen in the oocyte and during spermatogenesis, although there is also a small contribution of gonosomal mutations and mutations that arose in the early development of the child [44].We only select probands in which the parental origin has been determined for over 90% of germline mutation calls.The data set (see Table 2 below) includes cell types with varying cell division rates, ranging from postmitotic neurons to rapidly dividing epithelial cells from small bowel and colon [76].We excluded data from donors with cancers in the focal tissue, and past and current smokers from the lung data set [21].Download links to the data are provided in S1 Data.
We annotated each SBS substitution with its local context, i.e., the 5' and 3' flanking nucleotides, using the GRCh37 genome reference [77].Using this annotation, we classified each substitution into one of the 192 types (there are 4 3 = 64 different 3-mer contexts, and for each of them, there are 3 possible substitutions).This classification was then collapsed to give the standard 96 categories, which are strand-invariant [6,7].

Mutational signatures attribution
First, we describe how we decompose a set of mutations as a linear combination of 79 COS-MIC version 3.3 SBS mutational signatures [10].For any mutation of SBS type z, the probability that it is attributed to a given signature class s is given by Bayes' formula where P(z) = S s P(z|s)P(s).P(z|s) is known, but we need to infer the set of loadings P = {P(s)} from the data.The log-likelihood of the loadings P given the observed mutations fz j g N j¼1 is given by where P(z j ,s j ) is the joint probability of mutation of type z j and the underlying signature s j (a hidden variable).We maximize the likelihood iteratively using the expectation-maximization method (EM, [78]), under the constraint of normalization, S s P(s) = 1, and P(s)�0 for all s.Initial condition: We initialize the optimization with a uniform distribution.E step: Provided the distribution at iteration i, P i = {P i (s)}, we compute the pseudo-log-likelihood of the set of loadings P, M step: We find the distribution in the next iteration, P i+1 , by maximizing the pseudo-loglikelihood Q(P|P i ) under the constraint S s P(s) = 1, i.e., where Z is a Lagrange multiplier.In this way, we find that and the value of the multiplier can be found by applying the condition of normalization, Z ¼ P s P N j¼1 P i ðsjz j Þ.We apply regularization to avoid overfitting.Specifically, we impose sparsity on the distribution P i by introducing a Dirichlet prior over the set of distributions P, PðPÞ / e À b P s logPðsÞ ; ð21Þ where β�0 parametrizes the degree of sparsity.As discussed in the main text, the COSMIC catalog of signatures is based on cancer samples, and therefore some of the processes that shape the mutational spectra in noncancerous cells may be missing from this data set.It may lead to errors in signature attribution: mutations corresponding to undefined signatures will be incorrectly assigned to COSMIC signatures.One strategy to overcome this limitation might be to impose a high degree of sparsity of signature loadings, leading to a small number of attributed signatures and avoiding spurious signature attributions.However, because signature SBS5 comprises many mutation types, strong regularization could lead to the misassignment of many mutations to this signature, and consequently, its loading may be overestimated.Given our focus on signature SBS5, we instead chose to apply a weak regularization that avoids this problem, but at the potential cost of spuriously assigning mutations to other cancer-specific signatures that do not in fact operate in the germline.With this regularization approach, the distribution P i+1 in the next iteration takes a closed form.Taking advantage of the fact that a Dirichlet prior is conjugate to multinomial likelihood (16), it can be shown [79] that where the normalization factor reads Z i b ¼ P s maxf0; P N j¼1 P i ðsjxÞ À bg.This expression implies that if the estimated loading of a signature s is too low, it is set to 0.
In what follows, we extend this signature attribution method to infer age-dependent signatures in data consisting in mutations from donors of varying ages.First, we preform a linear regression on the number of mutations per haploid genome, y, on age, x y ¼ ax þ b: ð23Þ We find a and b for each cell type, using the scipy package [80] (see S4 Fig) .Next, we assign the mutational signatures that make up the slope and the intercept.Namely, we assume that independent distributions govern the signature compositions of the age-dependent and ageindependent mutations and denote them P a = {P a (s)} and P b = {P b (s)}, respectively.The expected proportion of mutations of type z at age x is then given by X s PðzjsÞðaP a ðsÞx þ bP b ðsÞÞ; ð24Þ where we sum over COSMIC SBS signatures.
In order to pose the problem of estimating the two components, P a and P b , we assign to 2 hidden variables to each mutation j: the signature to which the mutation is attributed, s j , and an indicator variable I j , which equals 1 for age-dependent mutations and 0 for age-independent mutations.In these terms, the log-likelihood of the data given the two sets of loadings is I j logðPðz j js j ÞP a ðs j ÞÞð1 À I j ÞlogðPðz j js j ÞP b ðs j ÞÞ: ð25Þ We estimate the loadings using expectation-maximization, as previously described.In this case, the pseudo-log-likelihood is where the two sets of membership probabilities P i a=b ðsjz j Þ are computed analogously to Eq 18. Inclusion of the age-dependent component leads to consistently higher likelihoods in all the data sets examined and lower Bayesian information criterion values in all data sets except for muscle and paternal mutations (S5 Fig), for which we find few differences between the agedependent and age-independent distributions.We assess the uncertainty of our estimates by bootstrapping.For each data set, we resample cells (or trios in the case of pedigree data) with replacement 1,000 times.For each replicate, we estimate the slope and intercept of the age dependence and infer the corresponding signature attribution.We test for the significant presence of a given signature by performing the inference with weak regularization (i.e., with β = 0.1).This choice of regularization strength does not affect the main components and only zeroes out the components of the decomposition with very low membership probabilities (see S6 Fig), as specified by Eq 22.We include a signature s from a given data set if the probability P a/b (s) is nonzero in at least 95% of bootstrap replicates; in the figures, we denote the residual contribution of signatures that do not meet this condition with an asterisk.

Fig 1 .
Fig 1. Mutagenic consequences of DNA damage.(A) Interplay of DNA damage and repair during the cell cycle.DNA damage leads to lesions at rate u per base pair.Lesions are repaired at rate r and lead to mismatches with probability �, due to the misincorporation of nucleotides by the DNA polymerase used in repair.Mismatches are resolved at rate q, resulting in the incorrect base pair and a mutation with probability 1-p.(B) Consequences of DNA replication.Replicating DNA over a lesion requires translesion synthesis.This process is not always accurate: it causes an error and a mutation in one of the 2 daughter cells with probability R (assuming that the lesion is repaired in the next cell cycle, i.e., that r�ϕ).Unresolved mismatches cause a mutation in one of the 2 daughter cells.(C) The predicted number of mutations, m, in a genome length l at age x contributed by the different mechanisms.The genome length is denoted by l and the rate of cell division by ϕ. https://doi.org/10.1371/journal.pbio.3002678.g001

Fig 2 .
Fig 2. Model of the accumulation of CpG transitions.(A) Model of the consequences of spontaneous deamination of methylated cytosines during the entire cell cycle[23].Double-stranded cytosines deaminate at a constant rate u d per base pair and the resulting mismatches are repaired at a rate r.With probability �, the mismatch resolution is incorrect and leads to a mutation.We assume that the cell divides immediately after DNA replication and treat the two processes as simultaneous (occurring at a rate ϕ).Unresolved mismatches at cell division lead to a mutation in one of the two daughter cells.We provide the prediction of the model for the number of mutations m at a given age x.Inefficiently repaired mismatches accumulate with the cell divisions, which occur at rate ϕ, and repair errors accumulate with absolute time, independent of cell divisions.The number of methylated cytosines is denoted by l. (B) An alternative source of CpG transitions is the deamination of methylated cytosines in single-stranded DNA (at a rate of u s ) during DNA replication.The model assumes that a deamination immediately preceding the polymerization of the second strand is not repaired and leads to a mutation in one of the daughter cells.The expected number of mutations is proportional to the rate of cell divisions and the time of transient single-strandedness Δt.

Fig 3 .
Fig 3. Clock-like mutational signatures in different cell types.(A) Age-dependent signature attribution of mutations in neurons.Shown is the increase in the number of mutations with age (reported per haploid genome); each point corresponds to a single donor.(B) The relative contributions of different signatures vary with age.The decomposition of the mutation spectrum into age-dependent (blue, C) and constant signature distributions (yellow, D).Mutational signatures are indicated by their COSMIC label; asterisks indicate unattributed signatures (see Methods).(E-L) The number of mutations assigned to clocklike signatures, SBS1 (red) and SBS5 (turquoise) in: (E) neurons, (F) maternal germline mutations, (G) paternal germline mutations, (H) smooth muscle from bladder, (I) liver hepatocytes, (J) lung epithelium, (K) small bowel epithelium, (L) colon epithelium (*for this data set, in which the intercept is negative, we assume both signatures increase with age).Throughout (E-L), shaded areas represent 95% confidence intervals, estimated by bootstrapping (see Methods).Underlying data for this figure can be found in S2 Data.https://doi.org/10.1371/journal.pbio.3002678.g003

Fig 4 .
Fig 4. Signature attribution for germline mutations.(A) Effect of maternal age on mutations assigned to the maternal germline in pedigree data.(B) Effect of paternal age on mutations assigned to the paternal germline in pedigree data.Paternal germline mutations (purple) accumulate at similar to rate to somatic mutations in testis seminiferous tubules (blue).(C-F) The decomposition of the mutation spectrum into age-dependent (C, E) and constant signatures (D, F) in maternal and paternal mutations; see (A) and (B) for the color code.The inset in (E) shows the decomposition of somatic mutations in the testis seminiferous tubules.SBS signatures are indicated by their COSMIC label; asterisks indicate the contribution of unattributed mutations (see Methods).(G-J) The distribution of SBS types reconstructed using the signature attributions (C-F).SBS types are grouped by the 6 substitution types (see J for color code) and sorted alphabetically by the sequence context (ACA to TTT).Underlying data for this figure can be found in S2 Data.https://doi.org/10.1371/journal.pbio.3002678.g004

Fig 5 .
Fig 5. (A) Kinetics of the interplay of damage and repair in the general case.The two-state model used for CpG transitions accumulation is obtained by taking the limit q!0.(B) Example solution of the model for the number of lesions (left), mismatches (center), and mutations (right).https://doi.org/10.1371/journal.pbio.3002678.g005 ÞlogðPðz j jsÞP a ðsÞÞ þ P i b ðsjz j ÞlogðPðz j jsÞP b ðsÞÞ:ð26Þ

Table 2 . Data sets used in the analysis of clock-like mutation accumulation.
s PðzjsÞP i ðsÞ: