Rate volatility and asymmetric segregation diversify mutation burden in mutator cells

Mutations that compromise mismatch repair (MMR) or DNA polymerase exonuclease domains produce mutator phenotypes capable of fueling cancer evolution. Tandem defects in these pathways dramatically increase mutation rate. Here, we model how mutator phenotypes expand genetic heterogeneity in budding yeast cells using a single-cell resolution approach that tallies all replication errors arising from individual divisions. The distribution of count data from cells lacking MMR and polymerase proofreading was broader than expected for a single rate, consistent with volatility of the mutator phenotype. The number of mismatches that segregated to the mother and daughter cells after the initial round of replication co-varied, suggesting that mutagenesis in each division is governed by a different underlying rate. The distribution of “fixed” mutation counts that cells inherit is further broadened by an unequal sharing of mutations due to semiconservative replication and Mendelian segregation. Modeling suggests that this asymmetric segregation may diversify mutation burden in mutator-driven tumors.


Introduction 26
All tumors contain genetically divergent cells spawned by the evolutionary processes of 27 mutation and selection. In some tumors, genetic heterogeneity arises from a "mutator 28 phenotype" 1 due to mismatch repair (MMR) defects 2 or heterozygous exonuclease domain 29 mutations (EDM) affecting the leading or lagging strand DNA polymerases (pol), Polε or Polδ 3-9 . 30 Since MMR corrects polymerase errors, when MMR and EDM mutations occur together they 31 produce a dramatic increase in the number of unrepaired polymerase errors. The resulting 32 tumors rapidly evolve and possess "ultra-hypermutated" genomes. Yet a full understanding of 33 the relative contributions of mutagenesis and selection to the rise of heterogeneity within these 34 tumors remains elusive, since cells with more mutations tend to adapt more readily. 35 A key unanswered question is whether the mutation rate is constant within mutator cell 36 populations. The two most common ways of measuring mutation rates are fluctuation analysis 10 37 and mutation accumulation lines 11 . Both assume a uniform mutation rate and report the average 38 of hundreds or thousands of cell divisions. However, in recent years, evidence has emerged 39 that mutagenic processes may vary from one division to the next. Kataegis and chromothripsis, 40 for instance, sharply increase mutation burden in a single cell division [12][13][14] . Indirect evidence for 41 highly mutagenic sub-populations of cells also comes from studies of yeast exposed to 6-42 hydroxylaminopurine or AID/APOBEC cytosine deaminase. Drug-resistant mutants in these 43 studies had substantially higher mutation burdens than non-selected isolates from the same 44 population 15 . More recently, limited single-cell propagation of human cancer cell lines coupled to 45 whole genome sequencing revealed broader than expected variation in mutation rate in closely 46 related subclones 16 . Observations such as these challenge the assumption that mutation rate is 47 constant and beg higher resolution studies of mutator cells. 48 The asymmetrically dividing budding yeast, Saccharomyces cerevisiae, is ideal for 49 studying mutator phenotypes with high resolution. It encodes many of the same DNA replication 50 and mismatch repair genes found in humans. Yeast "daughter" cells can be separated from their 51 larger "mother" cell at each division by micromanipulation and then moved to defined locations 52 on an agar plate, forming a "single cell lineage". Whole genome sequencing (WGS) of cultures 53 derived from these cells permits the number of new mutations that arose in the mother cell at 54 each division to be counted. Moreover, the small size of the genome (12 megabases) makes it 55 cost effective to score enough cell divisions to see whether the distribution of mutation counts 56 conforms to that expected from a single underlying mutation rate. 57 We previously pioneered this approach with haploid mutator mother cells deficient in 58 Polε proofreading and MMR (pol2-4 msh6Δ) 17 . A single underlying mutation rate could not 59 explain the distribution of mutation counts from 87 divisions. However, the distribution did fit a 60 model with two underlying mutation rates that differed by 10-fold (0.4 and 4 61 mutations/genome/division). This led to a hypothesis of "mutator volatility" in which cells 62 assumed one of two mutator states as they passed through the cell cycle 17 . But since we only 63 scored mutations retained by the mother, we could not exclude an alternative hypothesis: that 64 polymerase errors sporadically segregated asymmetrically between mother and daughter cells, 65 either as mismatches at the initial division or as permanent, "fixed" mutations following the next 66 round of synthesis. Here, to distinguish between these two hypotheses, we sought to score all 67 replication errors that arose in individual cell divisions using more extensive single cell lineages. 68 Examination of the distribution of the full replication error counts from individual divisions 69 provided a way to test the mutator volatility hypothesis apart from the confounding influence of 70 segregation. At the same time, sequencing complete lineages gave us the means to determine 71 whether replication errors segregate equally on their way to fixation. 72

73
To confidently score replication errors arising on all nascent DNA strands from each 74 division, we devised a scheme that ensured that all mutations were observed in at least two 75 members of a single cell lineage. After moving each daughter by micromanipulation from the 76 founding mother cell, we isolated a sublineage of three additional cells to help score the number 77 of errors segregated to that daughter. These cells included the first and second granddaughter 78 (born to the daughter cell) as well as the first great-granddaughter cell derived from the first 79 granddaughter (Fig.1a). Errors segregated to the daughter as mismatches in the first division 80 segregate as fixed mutations in the next division when the daughter produces the first 81 granddaughter. Mutations retained by the daughter after that segregation event will be inherited 82 by the second granddaughter, forming what we call the "Da" segregant group. Mutations 83 segregated to the first granddaughter will be inherited by the great-granddaughter, forming the 84 "Db" segregant group. In theory, the Da and Db segregant groups represent half of the errors 85 made by the mother cell during a given division. The remaining errors, retained initially by the 86 mother as mismatches, segregate between the mother and her next daughter as fixed 87 mutations in the next division. Fixed mutations segregated to that daughter will be uniquely 88 present in the next sublineage, forming the "Ma" segregant group. Mutations retained by the 89 mother will be found in all later sublineages, defining the "Mb" segregant group. After colony 90 formation and WGS, a full error count for a given division can be determined by simply summing 91 the number of fixed mutations in the Da, Db, Ma, and Mb segregant groups. With a complete set 92 of sublineages from the same mother cell, the full replication error counts from several 93 sequential cell divisions can be determined from the nested data (Extended Data Fig.1). By 94 requiring that all errors be observed in at least two members of the lineage, this approach 95 eliminates false positives due to sequencing errors or clonal sweeps within the cultures. 96 We initially began our experiments with the pol2-4 msh6Δ haploid strain used in the 97 previous study 17 . We found evidence for a more limited mutator volatility but were concerned 98 that lethality within some sublineages may have introduced a bias (see Supplementary  99 Information and Extended Data Fig. 2). To improve viability and the mutational signal, we 100 switched to using diploid yeast with a 10-fold higher mutation rate due to homozygous mutations 101 affecting Polδ proofreading and base-base mismatch repair (pol3-01/pol3-01 102 msh6Δ/msh6Δ) 18 lineage information allowed us to assign when the mutations arose using the logic described 112 above. In addition, we visually inspected the variant sites in all genomes from a given lineage 113 using the Integrative Genomics Viewer, which allowed us to detect discrepancies in the lineage 114 order or whether mutations had been incorrectly assigned (see Methods). We tallied the full 115 replication error counts from each division and determined whether the distribution could be 116 explained by a single underlying mutation rate. 117 Mutagenesis has been modeled for more than 70 years 18-20 with the Poisson distribution, 118 which is a discrete probability distribution of the number of expected independent events 119 occurring within a defined interval, assuming a constant rate (λ). A simple test of whether a 120 distribution matches a single Poisson is to calculate the index of dispersion (D ), which is equal 121 to the variance of the distribution divided by the mean (σ 2 /μ). The variance of Poisson 122 distributions always equals the mean, which results in a D of 1. The pol3-01/pol3-01 123 msh6Δ/msh6Δ mother cells committed an average of 276 (±37.7, standard deviation (σ)) 124 replication errors per division. This corresponds to a D of 5.15 (37.7 2 /276), which suggests that 125 the distribution does not conform to a single Poisson (Fig.1b). Two alternative explanations 126 failed to account for the overdispersion. For instance, we did not observe any relationship 127 between the mother's replicative age and the number of errors made by Polδ (Spearman's rank 128 correlation coefficient: 0.007209, p = 0.9604)(Extended Data Fig.5), nor did the number of 129 mutations correlate with the size of the scored genome, which differed between lineages due to 130 variation in sequencing depth and the number of members in each lineage (Spearman's rank 131 correlation coefficient: -0.0416, p = 0.7743)(Extended Data Fig.5). Instead, the broad 132 distribution of full replication error counts, free from the confounder of segregation, is consistent 133 with mutator volatility. 134 To better understand the nature of mutator volatility in pol3-01/pol3-01 msh6Δ/msh6Δ 135 cells, we used finite mixture modeling, which employs a maximum likelihood framework to 136 identify mixtures of two or more Poisson distributions that better fit the data. We also modeled 137 the data as a negative binomial (nb), which is a discrete distribution with separate rate (μ) and 138 shape parameters (θ) commonly used to interpret over dispersed count data. The rate 139 parameters λ and μ, for the Poisson and nb distributions, both define the mean number of 140 events. Since these models derive from different distributions, they cannot be directly compared 141 using standard statistical tests. Non-nested models such as these, however, can be evaluated 142 with Akaike Information Criteria (AIC), which uses maximum likelihood to estimate the loss of 143 information of each model relative to the observed distribution. To prevent overfitting, AIC 144 penalizes models with more parameters. Lower AIC values correspond to a more parsimonious 145 fit; however, interpreting differences in raw AIC values can be enigmatic. Thus, we transformed 146 the raw AIC values to "Akaike weighted values", which conveys their relative likelihood 147 (Fig.1b) 21,22 . We found that the negative binomial model was the most likely (relative likelihood 148 of 0.9999), followed by the two-Poisson-mixture model (2.2 x 10 -6 ), and the single Poisson (4.3 x 149 10 -28 ) (Fig.1b). Similar results were obtained using Bayesian Information Criteria (BIC), which 150 imposes stronger penalties for overfitting. Thus, mutator volatility in pol3-01/pol3-01 151 msh6Δ/msh6Δ cells is more complex than just two distinct mutator states. 152 The superiority of the negative binomial model suggests that the mutator phenotype may 153 vary continuously. This rationale derives from the ability to describe a negative binomial as a 154 gamma-Poisson distribution (Fig.2a). The gamma function is a continuous, rather than discrete, 155 distribution. Here, it takes the same shape parameter (θ) as the negative bionomial and serves 156 as a conjugate-prior to define variation in the rate parameter λ of a mixture of Poisson 157 distributions. The variation in λ that creates a negative binomial occurs between replication 158 events at the same site, or a collection of sites such as a chromosome or genome. Having 159 complete lineage information provided an opportunity to test whether λ varies at a chromosomal 160 or genome-wide level. The distributions of mismatches segregated to mother (Mm) or daughter 161 cells (Dm) across all divisions were the same and fit a negative binomial (Fig. 2b). If λ varied 162 widely during the replication of individual replicons (the units of DNA replication on a 163 chromosome), this could introduce asymmetry in the number of errors on sister chromatids, 164 which would then propagate to the daughter and mother cells (Fig. 2d). Consequently, Dm and 165 Mm from the same division would be free to vary within the observed negative binomial 166 distribution. Alternatively, if the genome-wide value for λ varies between cell divisions, a single 167 mutation rate would govern mismatch formation for both the mother and daughter genomes 168 (Fig. 2e). Dm and Mm would co-vary within the constraints of the corresponding Poisson 169 distribution. To distinguish between these two hypotheses, we first compared the correlation of 170 mismatches segregated to mother and daughter cells to simulated data generated under the 171 constraints of the two models. While no correlation was seen between Dm and Mm in the 172 simulated data from the first model (R 2 =0.001), similar correlations were observed for both the 173 simulated data from the second model (R 2 =0.47) and the actual data (R 2 =0.37). This 174 correspondence in the number of mismatches segregated to mother and daughter cells 175 extended down to the level of chromosomes (Fig. 2f). The R 2 values are lower than typically 176 seen with strong correlations, but as our modeling shows, this is expected since both X and Y 177 values are randomly drawn from a Poisson distribution. As a second test of the hypotheses, we 178 also performed 10,000 simulations of how each model would affect the distribution of full 179 replication error counts from 50 divisions (Fig. 2g). With the first model, the simulated index of 180 dispersion (3.28 ± 0.66, σ) was substantially less than observed with the actual data (D =5.15), 181 while the second model produced a good match (5.54 ± 1.12, σ). Together, these analyses 182 strongly suggest that the source of mutator volatility is variation in the genome-wide mutation 183 rate from one division to the next. To understand the potential implications of our findings for mutator-driven cancers, we 224 first focused on how the Poisson-binomial process would influence the heterogeneity of 225 mutation burden within a dividing population of tumor cells. Assuming a constant mutation rate 226 comparable to pol3-01/pol3-01 msh6Δ/msh6Δ yeast, the expected distribution of simulated 227 mutation counts in human cells after one division (D = 50) was far broader than in yeast (Fig.3f) 228 and persisted through 30 simulated divisions (Fig. 3g,h). Adding a comparable level of volatility 229 to the mutator phenotype further increased the simulated dispersion (D = 82) (Fig.3f). Using the 230 Poisson-binomial model, we simulated a range of mutator phenotypes observed in cancer cells 231 and found a linear relationship between mutation rate and predicted dispersion. For instance, 232 mutation accumulation in HCT116, the well-known MLH1 mutant colon cancer cell line, 233 increases from 48 to 190 mutations/haploid genome/division upon introduction of a 234 heterozygous POLE proofreading-deficient allele 9 . In these cells, the predicted index of 235 dispersion expanded from 3.4 to 10.8 (Fig.3i). Even greater heterogeneity may arise in human 236 cancers when more potent POLE mutator alleles occur in combination with MMR 237 deficiency 5,7,23,24 . Thus, the fundamental Poisson-binomial process of asymmetric segregation 238 has the potential to dramatically expand the diversity of mutation burdens present among a 239 population of human mutator cells. constant rate but segregate asymmetrically on the way to fixation. The design of our single cell 250 pedigrees ensured at least two independent biological observations for each mutation, which 251 allowed us to confidently assign more than 13,000 mutations to fifty divisions. From the resulting 252 mutation count data, we found strong evidence that both mutator volatility and asymmetric 253 segregation significantly expand genetic heterogeneity in pol3-01/pol3-01 msh6Δ/msh6Δ yeast. 254 Historically, mutagenesis has been modeled with the Poisson distribution, which 255 describes the probability of the number of independent events per unit time given a constant 256 rate. The observed distribution of full replication error counts of mutator cells, free from the 257 influence of segregation, best a fit a negative binomial and not a single Poisson (Fig.1b).

258
Negative binomials are equivalent to a continuous mixture of Poisson distributions whose rates 259 vary according to a gamma distribution (Fig.2a). This suggests that mutator volatility may create 260 a continuum of mutation rates rather than discrete mutator states. We explored the idea that 261 mutation rate varies from one division to the next by simulating the number of mismatches 262 segregated to mother and daughter cells (Fig 2d,e) and the dispersion of full replication error 263 counts expected from small cohorts of cells (Fig.2f). Both simulations closely matched the 264 observed data, consistent with the hypothesis that mutator volatility derives from continuous 265 variation in mutation rate between divisions. Mutator polymerases do not operate as a closed 266 system. They interface with a myriad of other replication components and metabolites, such as 267 dNTPs, that influence their fidelity 25,26 . Variation in the timing and duration of perturbations to 268 these interactions may produce a continuum of rates. The observed overall mutation rate that 269 cells exhibit represents a composite of mutation rates at all sites within the genome. 270 Conceivably, the change in replication fidelity could be localized to certain parts of the genome 271 in a given division. But if so, our data suggests, that the nascent strands from each pair of sister 272 chromatids in the affected region must be equally influenced by the change in rate (Fig.2c,f). 273 The asymmetric inheritance of mutations observed in mutator cells results from the 274 fundamental processes of semi-conservative replication and Mendelian inheritance acting in 275 concert. Current models of mutation accumulation generally ignore the potential for this synergy 276 to expand genetic heterogeneity, although there are exceptions. John Cairns proposed a far 277 more extreme asymmetric inheritance of mutations in the "Immortal Strand Hypothesis" in which 278 stem cells always segregated away newer DNA duplexes with fixed mutations 27 . In keeping with 279 this hypothesis, a recent computational analysis of human somatic variants argued that the high 280 variance of mutation burden in adult stem cells with age supports a preferential inheritance of 281 ancestral strands 28 . A second study from the field of evolutionary biology examined the potential 282 influence of disparate mutagenesis of leading and lagging strand synthesis to promote variable 283 evolutionary trajectories from the same cell population 29 . Our findings here demonstrate that, in 284 the context of a mutator phenotype, the normal process of semi-conservative replication and 285 Mendelian inheritance has the potential to create unequal sharing of mutations. For every cell 286 that inherits disproportionately more mutations there will be another cell with fewer mutations. 287 The predicted impact of this process on the variation in mutation burden is larger in human cells 288 than in yeast due to the vast differences in chromosome length, and the correspondingly larger 289 number of fixed mutations per chromosome. However, with longer chromosomes comes an 290 increased likelihood that sister chromatid exchanges (SCEs) may mitigate the asymmetry. 291 SCEs clearly to do not homogenize mutation burden in diploid mutator yeast cells as half of 292 cells either received all or none of the new fixed mutations for a given chromosome (Fig. 3c).

293
This finding is in keeping with recent evidence from a sensitive Next Generation should address both the prevalence of SCEs and the asymmetric inheritance of mutations. 302 Our simulation of a mutator-driven tumor rapidly generated substantial intra-tumoral 303 genetic heterogeneity during expansion (colored lines, Fig.3h) compared to a population in 304 which mutations accumulated by a simple Poisson process (black line, Fig.3h). The associated 305 variability in mutation load may be relevant to cancer evolution. Early during tumorigenesis the 306 subpopulation of cells that inherit disproportionately more mutations may adapt more readily. 307 With elevated mutation rates, polyclonal adaptation is almost certain. The unifying feature of 308 these adapted cells is a high mutation burden. As mutation burden mounts and mutator cells 309 contend with increasingly strong negative selection pressure due to immune surveillance and 310 negative epistatic interactions 32,33 , adapted cells that inherit fewer new mutations due to 311 asymmetric inheritance may be at a relative fitness advantage. Selectively increasing mutation 312 rate in mutator cancer cells could represent a novel therapy 25 . If, as a means of treatment, the 313 mutation rate of cancer cells is only transiently elevated to induce extinction, this subpopulation 314 may persist. Sustained elevation of mutation rate over many divisions of mutator cells may be 315 required to drive their extinction. 316

Single cell lineage isolation 333
To isolate pol2-4 msh6Δ lineages we dissected AH2801 tetrads on SC-Leu-Ura selective 334 media and chose one germinating spore per plate to serve as the founding mother cell. to allow mating. Upon isolation of a zygote, the first or second daughter was used as the 339 founding mother (M) for the lineage. Mothers were placed at an isolated location and we 340 separated daughter cells (designated Dn, Dn+1, etc.) from the mother as they were generated 341 and moved them to select areas 5 mm apart on the plate. We repeated the procedure to obtain 342 each daughter's first daughter (GD.1, Fig.1b), second daughter (GD.2), and first granddaughter 343 (GGD, born to GD.1). This strategy was repeated for each daughter up to either the 20 th division 344 or the end of the mother's replicative lifespan, whichever occurred first. In a typical experiment, 345 we pre-punched the agar with the dissecting needle at each drop-off location so that we would 346 always put the cell in a defined place, making it easy to later find the cell for inspection and 347 manipulations. We isolated lineages over the span of a week by performing rounds of 348 dissections every 90-120 minutes. Only a few cells on a plate were moved in any one round, 349 and then, only one cell at a time. We noted the timing of each round of bud dissections. We 350 incubated plates at 30°C between dissections. At the end of the day, plates were wrapped in 351 parafilm and stored overnight at 4°C. When plate dissections were concluded, we incubated 352 each plate an additional 48 hours at 30°C to allow colonies to fully develop. Prior to sequencing, 353 the pol3-01/pol3-01 msh6Δ/msh6Δ and pol2-4 msh6Δ genotypes were confirmed by previously 354 described allele-specific PCR assays 34 . 355

Genome sequencing 356
We performed whole genome sequencing of yeast cultures as described 34  variants should be present in 100% of reads. Setting the cut-off at 0.8 accommodates sites with 377 low read depth and one sequencing error. For pol3-01/pol3-01 msh6Δ/msh6Δ diploids, we used 378 a minimum read-depth of 18 for all strains and a variant frequency cut-off of 0.22. With a read 379 depth of 18, clonal heterozygous variants in diploid cells have a false negative rate of 6.1 x 10 -5 . 380 With 1000 mutations we have a 6% chance of having 1 false negative in a genome. We filtered 381 the above results to remove variants present in the parental strains as well as recurrent 382 sequencing artifacts. A small number of variants (<0.1%) could be reliably scored with the 383 above parameters but fell below a quality threshold for a subset of genomes. These were 384 manually curated for inclusion. We detected these by visually inspecting the BAM files for all 385 strains in a single cell lineage at the same time using the Integrated Genome Viewer (IGV). 386

Scoring of mutations and detection of assignment errors 387
We show up in all subsequent daughters (Dn+3, Dn+4, etc) and their offspring. Any deviation from 403 this pattern of inheritance indicates an "assignment error" has occurred, and that a cell was 404 inadvertently placed in the wrong position in the lineage. In the Supplementary Information we 405 describe two such cases. The divisions encompassing these strains were censored from the 406 analysis. Below we describe how these errors arise and are detected to illustrate the reliability of 407 the method. 408 One possible assignment error could occur at dissection when the daughter and mother 409 cells both divide before the next round of dissection. On the basis of size, the first daughter (Dn) 410 can be easily distinguished from the mother, the second daughter (Dn+1), and her own 411 daughter (GDn.1). Usually Dn+1 and GDn.1 can also be distinguished because Dn+1 buds 412 before GDn.1. However, in rare cases Dn+1 and GDn.1 are adjacent and similarly sized. If 413 Dn+1 is moved in place of GDn.1, we will have a sublineage consisting of Dn, Dn+1, GDn.2, 414 and GDn+1.1 (instead of Dn, GDn.1, GDn.2, and GGDn.1). Every sublineage should normally 415 contain subsets of mutations from different divisions (Da and Db mutations from the "n" division; 416 Ma mutations from the "n-1" division; and Mb mutations from the "n-2" division). In this 417 sublineage, the Ma segregant group mutation count will be 0, since there are no new mutations 418 that will be shared by these four colonies. However, a substantial subset of the mutations 419 assigned to the Db segregant group will also be found in later sublineages indicating that they 420 are not Db mutations but Mb mutations from a later division. The other half of what appear to 421 be Db mutations will in fact be Ma mutations from a different division. Added confirmation of the 422 dissection error comes from the analysis of the next sublineage, which will consist of GDn.1 (not 423 Dn+1 as it should be), GGDn.1, GGDn.2, GGGDn.1 (great-great-great granddaughter 1). There 424 will be 0 Mb mutations in this sublineage since all of these cells are directly descended from Dn. 425 These problematic cell divisions would be censored because we lack key lineage members 426 necessary to obtain a full replication error count. inheritance to flag this as an error. We don't think this is a common problem given the 442 correspondence between mismatches segregated to the mother (Mm) and daughter (Dm) cells 443 illustrated in Fig.2d Division 8, is also associated with a Ma/Mb pair with high mutation counts (120,65), leading to 454 the conclusion that this division had a high mutation rate. 455

Statistical modeling 456
We grouped the fixed mutation counts from the above branch points into Da, Db, Ma, 457 and Mb segregant groups to determine their distributions. We also joined all segregant groups 458 into one larger group to examine the distribution of fixed mutation counts across all cell 459 divisions. To determine the distributions of mismatches segregated to the daughter (Dm) and 460 mother (Mm) cells, we first summed the Da and Db or Ma and Mb fixed mutation counts from 461 each division. We also combined these two sets into one group to view the distribution of 462 mismatches across all cell divisions. To determine the distribution of total polymerase errors per 463 division, we summed all fixed mutations from individual divisions (Da+Db+Ma+Mb). We 464 considered two common approaches for modeling over dispersed count data: the Poisson 465 mixture distribution and the negative binomial distribution. From this formulation, we see that the full density of the distribution is decomposed as a sum of 474 the scaled Poisson densities. In (1), represents the prior probability that a given count 475 measurement will be generated from the kth Poisson component distribution, parameterized by 476 λ . Since a given count measurement could have been generated from any of these K 477 components, we average over their densities based on their prior probabilities to get the full 478 density of that count. 479 The negative binomial distribution can be specified by the following probability mass where μ is the rate parameter and θ is the shape or dispersion parameter. As θ tends towards 483 zero, the variance increases. As θ → ∞, the negative binomial reduces to a Poisson 484 distribution. 485 We implemented these principles using a single R script (FMM.R, see Github link 486 below).To fit Poisson mixture models we used the flexmix R package in R v3.5.3 41 . To fit 487 negative binomial models we used the glm.nb function of the MASS R package 42 . Goodness of 488 fit testing of the models was performed using both Akaike information criterion (AIC) and 489 Bayesian information Criterion (BIC) in R. Although these two approaches score fit in slightly 490 different ways, BIC returned results consistent with AIC and we thus report only the more 491 commonly used AIC scores. We scored each tested distribution against up to 4 parameters. We 492 reported only up to the number of parameters that improved model fit. Lower raw AIC values 493 indicate better fit; however, the relative differences are not immediately intuitive and so we 494 calculated Akaike weighted values as described 21,22 . To illustrate this approach, the AIC values 495 in Fig.1b were 637, 537, and 511. The first step in getting weighted AIC values is to determine 496 ΔiAIC: the difference between each AIC value and the AIC with the lowest value (so for these 497 numbers: 126, 26, 0). The likelihood of each is then calculated by exp(-1/2 x ΔiAIC). The 498 weighted AIC value for a given model is its likelihood divided by the sum of all competing 499 likelihoods. From these calculations the weighted AIC values are 4.3e-28 (P, k=1), 2.2e-6 (P, 500 k=2), and 0.9999978 (nb), respectively. Thus, the negative binomial model is far more likely 501 than the other two models to account for the observed data. Mixture model graphs were 502 constructed using the ggplot2 package R 43 . Spearman rank correlation coefficients were 503 calculated using the Scipy Stats package in Python and graphs generated with Seaborn 0.9. 504

Simulation of negative binomial models 505
We wrote a Python script (Fig2a-d.py, see Github link below) to simulate the expected 506 correlation between Dm and Mm under two distinct models of mutagenesis (Fig.2). The script 507 uses the θ (60.42) and μ (138) parameters estimated by glm.nb for the negative binomial model 508 of mismatches segregated to mother (Mm) or daughter (Dm) cells (see FMM.R). (Note that 509 glm.nb actually returns the natural log value for μ (in this case 4.927), which must be 510 exponentiated (e 4.927 ) to get 138). In the first model, we assumed that the negative binomial 511 distribution was created by variation in mutation rate along chromatid pairs, so that upon 512 segregation, Dm and Mm from the same division were free to vary within the predicted negative 513 binomial distribution. To simulate this process with Scipy.stats.nbinom.rvs, we converted the θ 514 and μ shape parameters to the n and p inputs (see script for details) for nbinom.rvs and then, for 515 each division, we selected two random values from the distribution to represent the Dm and Mm 516 counts. In the second model, we assumed that the negative binomial was created by a gamma 517 distribution of λ values for a series of Poisson processes acting in different cell divisions. We 518 used Scipy.stats.gamma.rvs to simulate λ values from a gamma distribution with shape and 519 scale parameters derived from those of the negative binomial. The shape parameter for the 520 gamma distribution is simply equal to θ. With variance (v) equal to μ 2 /θ, the scale parameter is 521 equal to v/μ. With a random λ from the gamma distribution as an input for 522 Scipy.stats.poisson.rvs, we selected two values from the associated Poisson distribution to 523 serve as Dm and Mm counts for each division. To examine the relationship between Dm and 524 Mm in these different models and the actual data, we performed linear regression with 525 Scipy.stats.linregress and visualized the data and regression line using Seaborn 0.9 regplot. 526

Simulation of gamma-Poisson-binomial process 527
We and then divided these values by the total length of unmasked DNA in the haploid genome. The 532 rate of mismatches per haploid genome (69 mismatches/haploid genome/division for pol3-533 01/pol3-01 msh6Δ/msh6Δ cells) was then multiplied in each case by these fractions to obtain 534 the per chromosome rate of mismatch formation. These values were used as input for 535 scipy.stats.poisson.rvs to simulate the number of errors per chromosome in a single division. 536 We created two independent entries per chromosome to model the diploid genome. To mimic 537 the binomial process of Mendelian segregation, we then multiplied the number of simulated 538 errors on each chromosome by a randomly chosen 1 or 0. Finally, we summed the mutation 539 counts from all chromosomes to obtain the total number of fixed mutations per cell. To create a 540 gamma-Poisson-binomial model, we selected a value for lambda at each division from the 541 gamma distribution described in Fig.2 rather than using a constant rate for mismatch formation. 542 As a control we performed the above simulation without the binomial process, using the rate of 543 fixed mutations per haploid genome (34.5 fixed mutations/haploid genome/division). We used 544 the same approach for the human simulations except that we multiplied the fraction of each 545 human chromosome of the total genome (GRCh38) by a mismatch rate comparable to that 546 observed with pol3-01/pol3-01 msh6Δ/msh6Δ yeast: 69 mismatches/ haploid yeast 547 genome/division x (3.03 x 10 9 bp/human haploid genome / 11 x 10 7 bp/yeast haploid genome) = 548 1900 mismatches/human haploid genome/division. We compared the resulting distribution to 549 that from a Poisson distribution with λ equal to 950 fixed mutations/haploid genome. To simulate 550 the diversity in mutation burdens that this process generates, we summed the simulated 551 mutation counts for individual lines from 30 divisions. 552

553
Sequence data used to generate the findings of this study have been deposited in the 554 NCBI Sequence Read Archive (SRA), BioProject accession: PRJNA586886. 555 Scripts used to generate figures and perform statistical tests have been deposited to github: 556 https://github.com/idowsett/Asymmetric-segregation-of-polymerase-errors-and-rate-volatility-557 diversify-mutation-burden 558