Worm strains and maintenance
Bristol N2 was maintained as described [79] using E. coli HB101 as a food source. Worms were staged by picking at the L4 stage as identified by visual inspection of vulva shape. For larval germ-cell counts, young adults were transferred to fresh plates every 2 h for 8 h to produce several synchronized egg populations. The embryos were incubated for 21 h from the initial collection point and the larvae were dissected at approximately 2-h intervals, so that larvae used were collected between 21 and 36 h after being laid. Germ cells were identified by staining for PGL-1 [80]. In the course of the same experiment, populations were set aside and were not sacrificed for germ-cells counts but were observed at 2 h intervals on the day they were expected to reach L4; it took 54 h from the time of egg laying for 90 % of the population to have reached mid-L4 (based on scoring n ⩾ 50 worms at each time point).
Staining and imaging
For EdU pulse-chase experiments, worms were fed E. coli MG1693 that had been grown in minimal medium supplemented with glucose [81] and 75 mM of the thymidine analog EdU (C10337, Life Technologies, Grand Island, NY). Immediately following seeding, plates were stored at 4 °C. Plates were warmed to 20 °C prior to use. Worms were kept for either 15 or 30 minutes on EdU-labeled bacteria in the dark, returned to non-labeled bacteria in the dark for the period of the chase, and were fixed and processed as described [26] using 0.1 μg/ml DAPI to label DNA and 1:200 anti-PH3 antibody (9706, Cell Signaling, Beverly, MA) followed by Alexa 594-conjugated anti-mouse antibody (A21203, Life Technologies, Grand Island, NY) to label M-phase cells.
CYE-1 and PGL-1 stainings were performed by freeze-cracking dissected gonads or whole larvae, dehydration in acetone, 5-minute fixation in 4 % PFA, incubation with anti-CYE-1 antibody (a gift from Edward Kipreos) at 1:5 dilution or rabbit anti-PGL-1 antibody (a gift from Susan Strome) at a 1:500 dilution, and incubation with DAPI and Alexa 594-conjugated anti-mouse or anti-rabbit secondary antibodies.
All samples were imaged at ~0.3-μm z intervals with LSM 710 or 780 confocal microscopes (Carl Zeiss MicroImaging, Oberkochen, Germany), using a 63× objective.
Computational simulations
We developed a computational model of germ-cell cycling and differentiation as sperm or oocyte, and implemented it in C++. The same computational core is used for simulations of mutation accumulation and for fitting of cell cycle parameters to experimental EdU pulse-chase data.
Rules for cell movement and differentiation
In the MZ, a two-dimensional lattice is considered that has a long axis (corresponding to the distal–proximal axis of the gonadal arm) and a short axis that wraps around itself to form a hollow cylinder mimicking the shape of the gonadal arm. Only one cell can occupy a lattice point at any given time. When a cell divides, one daughter remains at the same location and one daughter needs to find a new position. If an empty lattice point exists in the same row the division occurred, cells in the row are pushed across the short axis so that the nearest empty point in the row is filled. Otherwise, if the next cell row has an empty position, the daughter cell is pushed forward to that row, and cells within the new row are displaced as necessary so that the empty position is occupied. If both the row in which the division occurred and the next row are full, the daughter is either pushed forward to the next row or sideways in the same row with equal probability and thereby displaces another cell. The same movement rules are then iteratively applied to this displaced cell and other cells that are subsequently displaced, until either an empty point is filled in the MZ or a cell is pushed out of the MZ. The randomness in simulated cell movement is inspired from the randomness observed in the orientation of cell division planes [27].
Within the meiotic zone (MeZ), eventual cell fate is either spermatogenesis, oogenesis, or apoptosis. Cells are drawn upon for these fates depending on their time of entry into the MeZ (cells are pushed out in first-in-first-out fashion), in a way that allows worms to maintain homeostasis of germ-cell numbers and to meet the experimentally defined development rate and reproductive schedule (as detailed below).
In some simulations (8, 9, 11–13 in Table 1 and Additional file 1: Table S1), the shape of the MZ was allowed to change with worm age to match experimental behavior. For rows that see their capacity diminish, thus requiring cell rearrangement, the same movement rules as above are applied.
Rules for timing of cell division
Cell cycle length is defined using a piecewise-linear function of position along the distal–proximal axis and of developmental time. The number of control points was kept to its smallest useful value to ensure computational tractability and avoid overfitting. For fits to experimental data, there was no temporal freedom given the short length of the EdU chase, and we used two or three spatial control points. For simulations of mutation accumulation, the number of control points we used is reported in Additional file 1: Table S1.
At each cell birth during the simulation, a time of next cell division is computed independently for the two daughters by sampling from a uniform distribution whose mean is determined by the piecewise-linear function described above and whose width is 1 % of cell cycle length. For simulations considering cells that can stop cycling within the MZ (see section “Rules for leaving the cell cycle”), the time at which the cell will enter G2 is computed using experimentally determined cell cycle phase indices (Fig. 3b). Both these times are entered into a priority queue that keeps track of the next event to take place in the simulation. If the cell is pushed forward before it has divided, the time to next division is scaled using the ratio of cycle lengths between the new row and the old row.
Time in the simulation moves forward by retrieving the next simulation event from the priority queue each time the previous event – a cell division and ensuing displacement events or a cell leaving the cycle – has been processed.
Rules for leaving the cell cycle
A fraction of cells within the PMZ leave the mitotic cell cycle but do not immediately proceed with meiosis; these cells have been referred to as pre-meiotic [27, 28, 45, 82]. The drop in M phase along the distal–proximal axis (Fig. 3b) is consistent with pre-meiotic entry in the proximal region (as previously reported [28]), as is the concomitant rise we observe in G2 DNA content.
To model the process by which cells leave the mitotic cell cycle while still in the MZ, we assumed that cells reaching G2 in the PMZ could make a decision to proceed with another mitotic cycle or to arrest in a pre-meiotic state. In simulations that took pre-meiosis into account, each time a cell reached G2 in the PMZ, the program checked whether the simulated local mitotic index was higher than the experimentally derived index; if it was, the cell was arrested at G2 and directed to a meiotic fate. Given the difficulty in ascertaining which particular cells are in pre-meiosis and which are not, we further assumed that cell cycle indices remained constant throughout the PMZ for the proliferative fraction. Resulting fits show an excellent match to experimental data (Additional file 2: Figure S4).
Our fits to experimental data taking pre-meiosis into account result in ~227 actively cycling cells within the population of ~257 cells comprising the MZ, with thus 30 pre-meiotic cells arrested in G2; note that at any given time a substantial fraction of the 227 cycling cells, found predominantly in the PMZ, will actually not get a chance to undergo another round of mitosis and could therefore also be considered pre-meiotic. Based on our best-fit simulations, the rate at which cells are pushed out of the MZ is ~20 cells per hour, in close agreement with an experimental measurement of that rate [45].
Initial conditions and result collection
For simulations of mutation accumulation, gonadal arms are seeded with a single progenitor cell. Each cell keeps a record of the number of divisions that link it to the progenitor cell, i.e. its pedigree depth. Depending on the purpose of the simulation, average pedigree depth is computed either from all cells leaving the MZ or from gametes that led to progeny production. In relevant cases, the pedigree-depth metric is adjusted to match variations in underlying biological assumptions. Specifically: 1) in immortal strand simulations, the pedigree depth of daughter cells that stay in the distal-most row is not incremented while the pedigree depth of daughter cells that are pushed forward is incremented by two (in such simulations, one daughter cell is always pushed forward upon cell division in the first row) and 2) when the mutation rate is assumed to depend on cell cycle length, the pedigree depth is incremented by 1/Ɣ, where Ɣ is the cell cycle length of the cell divided by the minimum cell cycle length over all the MZ (this normalization is applied so that pedigree-depth results are non-dimensional with respect to absolute cell cycle length).
For fitting to experimental cell cycle data, gonadal arms are seeded with a population of cells whose initial age within the cell cycle is taken from an exponential distribution, and whose overall cell phase distribution matched experimentally determined cell cycle phase indices (Fig. 3b). Simulations were pre-run for a period of 2 h in simulation time, at which points cells in S phase were marked as EdU-positive. Simulations were then further run for various amounts of time, and the distribution of cell cycle progression recorded at relevant simulated chase times for comparison with experimental data.
Computation of population growth rate
We also used our simulations of germ-cell cycling and differentiation to quantify the impact of slow germ-cell cycling on the overall population growth rate. For a given profile of cell cycle lengths along the MZ, we recorded the times at which cells destined to become oocytes were pushed out of the MeZ. Since fertilization occurs concomitantly with ovulation, this defined progeny birth times. We computed an average reproductive schedule based on 450 simulation runs, and used that schedule to define a transition matrix whose dominant eigenvalue yielded average population growth rate [83].
Derivation of experimental numbers required for optimizations
Optimizations require the input of numerical values for parameters that are not optimized and are thus derived from experimental data. Constraints used to enforce timely development and reproduction must also be derived from experimental data. We detail below how we used experimental data to set up optimizations, and how we applied constraints to the simulations.
Apoptosis
When considering gamete production, pedigree-depth optimizations need to take into account the fact that not all germ cells produced become gametes: a number undergo apoptosis (or contribute to growth of the MeZ without leaving it by the time sperm depletion stops reproduction). It has been reported that apoptosis is initiated as hermaphrodites transition to adulthood and switch to oogenesis [42], but the rate at which apoptosis occurs from that stage has not been fully determined: it is only known to be 50 % or more [42]. We extended this result and fully defined the apoptosis probabilities as a function of worm age. Because the rate of apoptosis is difficult to measure directly (counts of cells undergoing apoptosis at a given point in time do not readily translate to apoptosis rates), the idea we followed was to use the difference between distal cell influx into the MeZ from the MZ (inferred from cell cycle speed measurements) and proximal oocyte efflux (inferred from the reproductive schedule): after accounting for changes in MeZ size through developmental time, this difference provides the rate at which cells are eliminated. We implemented this idea using our simulations, relying on the following experimental data: measured cell cycle rates, MZ geometry (Additional file 2: Table S9), MeZ size at L4 + 1 and L4 + 3 (determined as 749 cells and 1077 cells, respectively, by subtracting MZ size from total germ-cell counts performed on n = 19 gonadal arms). We determined which apoptosis profile made it possible to match the experimentally characterized reproductive profile (Table 3), also adjusting the size of the simulated L4 MeZ (which is not completely filled at that stage) to 500 cells, which allows the first oocyte to be pushed out of the MeZ, thus initiating reproduction, at the appropriate time. The resulting apoptosis profile is shown in Table 2. In all subsequent simulations incorporating apoptosis, cells undergoing apoptosis were chosen stochastically, with a probability following this profile.
The effect of apoptosis on pedigree depth led us to define its temporal profile in the simulation as a function of total number of germ-cell divisions rather than elapsed time, for reasons detailed in the next two paragraphs. First, we note that the presence of apoptosis leads to an overall increase in gamete pedigree depth: for a given number of gametes to be produced, more germ cells need to be produced by the MZ if a number of these cells are fated for apoptosis instead of gametogenesis, which requires more cycling and thus a pedigree-depth increase. Second, because our experimental data show that cells leaving the MZ late in life are more likely to undergo apoptosis than cells leaving the MZ earlier (Table 2), taking apoptosis into account preferentially increases pedigree depth of late-produced gametes.
Third, we note that cell cycling in the MZ sets the pace at which reproduction proceeds in our simulations: faster germ-cell cycling leads to faster filling of the MeZ, faster pushing out of proximal MeZ cells once the MeZ is full, and thus faster oocyte maturation and reproduction. Therefore, if the apoptosis probability were defined as an increasing function of time, this could cause pedigree-depth optimizations to artificially favor a high speed of MZ cycling in development and early adulthood to compress the reproductive schedule to early adulthood; this would cause an unrealistically high rate of early reproduction. To avoid this behavior, we defined the apoptosis profile within the simulation as a function of total number of germ cells produced rather than as a function of developmental time. Put simply, before the first n
1 cell divisions have occurred, cells leaving the MeZ have an apoptosis probability of p
1; up to the next n
2 divisions, these cells have an apoptosis probability of p
2, etc. (where the p
i
are derived from experimental data as explained above). With this scheme, an overall speedup in the cell cycle leaves pedigree depth unaffected, which avoids an artificial pressure for cycling at high speed in young adults. For consistency, other simulation parameters that are dependent on time (MZ and MeZ geometry) and the position of temporal control points were also defined in terms of total germ-cell divisions that had occurred up to that point; only production constraints (detailed below) were defined as a function of elapsed time.
Timing of developmental stages
The temporal points that we used in experiments (L4, L4 + 1 day, and L4 + 3 days) were defined as developmental stages (worms were picked at mid-L4 based on vulva morphology, and used immediately or after 1 or 3 days; for simplicity mid-L4 is referred to as simply L4 throughout). For the purposes of our simulations, these stages needed to be expressed as total number of germ-cell divisions that had occurred in a gonadal arm (as discussed above) and as amount of elapsed time (to define production constraints). We estimated the number of cell divisions that had occurred by these stages using simulations of cell cycling in which the spatiotemporal cell cycle profile was set to its experimentally determined value (the number of cell divisions cannot be derived directly from cell counts because of cell loss to apoptosis and gametogenesis). We measured the amount of time that elapses from the time of egg laying to L4 as 54 h (see section “Worm strains and maintenance”). Given that the first sign of germ-cell proliferation is seen at 24 h, 30 h elapse between the onset of germ-cell proliferation and mid-L4. Overall, we defined the stages as follows: onset of germ-line proliferation: at 0 divisions or 0 h; mid-L4 (referred to as L4) at 400 divisions or 30 h, mid-L4 + 1 at 1200 divisions or 54 h, and mid-L4 + 3 at 2400 divisions or 102 h.
Production constraints
Pedigree-depth minimization favors a steep cell cycle gradient across the MZ, as shown by optimizations 1 and 5. A steep gradient can be achieved by fast cycling of proximal cells or slow cycling of distal cells. Physical limits on cell cycle speed must derive in part from the rate at which nutrients are processed and delivered to germ cells, and also from the speed at which the cells can replicate DNA and cellular structures; these limits were inferred from experimental data. We set the maximal cell cycle speed to the highest speed observed either during early development (3.4 h for L2–L3), or at any time of development and adulthood (2.8 h from L4 onwards). Slow cycling of distal cells results in a lower rate of cell production by the MZ. Simulations either incorporated constraints on the number of gametes produced by specific stages (following Table 3) and the total number of germ-cell divisions that had occurred by 30 h, 54 h, and 102 h (see above), which effectively placed a lower bound on the speed at which distal cells could cycle, or explicitly incorporated the cost of delayed reproduction caused by slow germ-cell cycling in the objective function being optimized. Overall, production constraints were such that the slowest possible cell cycle speed for distal cells was ~6 h at L4 (optimizations 10 or 11) and 16 h (optimization 11), or 32 h (optimization 10) at L4 + 1 day.
Geometry
For simulations that used an MZ geometry modeled directly after experimental data, we used measurements of numbers of cells per row at L4, L4 + 1 day, and L4 + 3 days. That geometry was linearly interpolated in time, and assumed to be constant before L4 and after L4 + 3 days. Details of parameters that change on a row-by-row basis and are linearly interpolated in time are given in Additional file 2: Table S9.
Image analysis
To quantify DNA, EdU, and CYE-1 contents cell by cell in intact gonads, we acquired three-dimensional confocal stacks of dissected tissue at short z intervals and used custom software to segment cells, i.e. to partition image pixels into distinct subsets that correspond to given cells (the software will be reported elsewhere and has been released as open source at [84]; image datasets can be downloaded [85]). Randomly chosen sample segmentations are shown in Additional file 2: Figure S5. To quantify the position of each segmented cell, we computed the geodesic distance to the distal end along the distal–proximal axis (using a principal curve computed as described [86]). To avoid artifactual attenuation of fluorescence intensity in cells deep in the tissue due to scattering and absorption, we only kept cells that had a direct line of sight to the microscope objective and thus exhibited minimal attenuation; such cells were identified by a metric we call top layer, defined as the relative cross-sectional area of their segmentation mask that projected to the top slice unhindered by masks of neighboring cells. An alternative method would have been to select stack top cells based on z position within the top nth percentile of z positions, where n can be adjusted stack by stack so that each stack contributes a given number of cells. Because attenuation is stronger when light travels through tissue than when it travels through an immersion medium, and because z variation throughout MZs was overall small (6 μm between MZ rows 1 and 10, n = 18), the top layer metric provided more accurate fluorescence quantification than the stack top metric (Additional file 2: Figure S6I, J; z position percentile adjusted so that both metrics selected the same number of cells per MZ, to ensure a fair comparison). DNA and EdU content were computed by summing all pixels within the cell, while nuclear CYE-1 contents were computed by summing pixels in a 0.4 × 0.4 × 1 μm box centered on the nucleus; DNA and CYE-1 contents were normalized so that the 10 % and 85 % quantiles mapped to 2C and 4C, respectively (DNA) or to 1 and 2 arbitrary units, respectively (CYE-1). Then 10 % and 85 % quantiles were chosen so that G1 and G2 peaks in the EdU-negative fraction of EdU pulse-fix experiments were correctly positioned at 2C and 4C. We further validated this normalization scheme on M-phase DNA contents across the full range of chase times (Additional file 2: Figure S6a–h). We also verified that the variation in CYE-1 signal along the distal–proximal axis was not an artifact of the deeper position of distal cells in image stacks (Additional file 2: Figure S6M).
Cell cycle fits
The aim of the cell cycle fits is to find the spatial cell cycle length profile that best fits experimental data. The overall procedure is to perform cell cycle simulations as described above, sampling free cell cycle parameters from a grid, and to report the set of parameters that provides the best fit to experimental data as measured using the DEMD or FLM metrics defined below. The free parameters are used to define total cell cycle length at one control point located at the distal end, one located at the proximal end of the MMZ, and, for simulations that encompass the PMZ, one located at the proximal end of the MZ. The cell cycle profile at positions in-between control points is linearly interpolated, as for pedigree-depth simulations. Unlike total cell cycle length, the relative lengths of G1, S, G2, and M can be directly computed from experimental EdU pulse data (0 h chase); assuming an exponentially decreasing cell age distribution f such that f(0) = 2 f(1) [87], we computed the relative length of G1 by solving F(x) = p
G1, where F(x) = 2 – 21-x is the fraction of cells younger than x according to the exponential age distribution, and p
G1 is the observed G1 phase index (and so forth for subsequent phases). These experimentally determined relative lengths are fixed in the simulation (numerical values are shown in Additional file 2: Table S9), and used in addition to total cell cycle length, which is defined by free parameters, to track cell progression through the phases of the cycle. Therefore, the free parameters that define total cell cycle length at their respective control points (for a total of two or three control points depending on the kind of simulation) fully define cell cycle behavior in the simulation. For each simulation run, a record is output that contains the value of the free parameters, and for each simulated cell its cycle phase, its DNA content (as computed from its simulated progression through S phase), its EdU content, and its position along the distal–proximal axis. This set of records is compared to experimental data using two independent metrics.
The two independent metrics that we used are DEMD, a new metric we developed that has the advantage of using all cells in the samples – thus providing information about all cell cycle phases and decreasing uncertainty in cell cycle parameter estimates – and FLM, a well-established technique that only makes use of the relatively small number of M-phase cells [46]. Briefly, DEMD measures the similarity between experimental and simulated DNA content histograms of EdU-positive and EdU-negative populations, while the FLM distance measures similarity between experimental and simulated FLM matrices. Algorithmic details for the FLM and DEMD metrics are given below.
Fitting using DEMD
Consider a series of EdU pulse-chase experiments across T different chase times. Suppose we quantify DNA content, EdU content, and spatial compartment C for each individual germ cell in our EdU pulse-chase experiments. It is then straightforward to generate a set of T × C × 2 DNA content histograms, where cells are partitioned based on chase time T, spatial position C, and EdU content (labeled or unlabeled). Define DEMD histograms as the set of histograms
$$ g = \left\{{g}_1, \dots,\ {g}_{T\times C\times 2}\right\},\ h = \left\{{h}_1, \dots,\ {h}_{T\times C\times 2}\right\} $$
Consider two sets of DEMD histograms g and h. Define the DEMD distance d
DEMD
between g and h
$$ {d}_{DEMD}\left(g,\ h\right) = {\varSigma}_in\left({g}_i\right)\ n\left({h}_i\right)\ {d}_{CEMD}\left({g}_i,{h}_i\right) $$
where n(.) gives the number of cells in a histogram and d
CEMD
is the circular earth mover’s distance [88]. Now, suppose g is drawn from experimental data and h(v) is drawn from simulations with cell cycle profile v. The goal of DEMD-based cell cycle fits is to perform the following minimization:
$$ {v}_{DEMD} = arg\ {min}_v{d}_{DEMD}\left(g,\ h(v)\right) $$
We found v
DEMD
via a grid search implemented in MATLAB.
Fitting using FLM
Consider a series of EdU pulse-chase experiments across T different chase times. Suppose we quantify cell phase, EdU content, and spatial compartment C for each individual germ cell in our EdU pulse-chase experiments. It is then straightforward to generate a T × C matrix that records the percentage of M-phase cells at chase time T and spatial position C that are EdU-positive. Define this T × C matrix as the FLM matrix p:
$$ p = \left\{{e}_{t,\ c}\right\} $$
Consider two FLM matrices g and h. Define the FLM distance d
FLM
between g and h:
$$ {d}_{FLM}\left(g,\ h\right) = {\varSigma}_{t,\ c}n\left({g}_{t,\ c}\right)\ n\left({h}_{t,\ c}\right)\ {\left({g}_{t,c} - {h}_{t,c}\right)}^2 $$
where n(.) gives the total number of M-phase cells used to compute the percentage. Now, suppose that g is drawn from experimental data and h(v) is drawn from simulations with cell cycle profile v. The goal of FLM-based cell cycle fits is to perform the following minimization:
$$ {v}_{FLM} = arg\ {min}_v{d}_{FLM}\left(g,\ h(v)\right) $$
We found v
FLM
via a grid search implemented in MATLAB.
Confidence intervals
We computed confidence intervals on v
DEMD
and v
FLM
via bootstrapping [89]. We performed bootstrapping in a way that each sample maintained the same number of gonadal arms at each chase time. Specifically, suppose we use an experimental dataset z composed of N gonadal arms. Suppose z is partitioned into T subgroups based on chase time:
$$ \begin{array}{l}z = \left\{{z}_1,\ {z}_2, \dots,\ {z}_T\right\}\\ {}{z}_i = \left\{{g}_1,\ {g}_2, \dots,\ {g}_{n(i)}\right\},\end{array} $$
where g
i
is a gonadal arm and where n(1) + n(2) + … + n(T) = N.
A bootstrap distribution for v
EMD
and v
FLM
was derived by resampling each z
i
independently and rerunning the grid-search minimization.
Pedigree-depth optimization
Our simulations of mutation accumulation are by nature stochastic, because of the randomness in cell movement and in cell cycling. Asking what cell cycle profile minimizes mutation accumulation thus requires minimizing a stochastic objective function, and deriving a range of parameters that perform reasonably well around that minimum. Compounding the difficulty of the problem, many of our simulations are performed under constraints on speed of development and reproduction that are subject to the same stochastic fluctuations.
We took a two-step approach to identify parameter sets that minimize mutation accumulation. First, we performed a grid search for parameters that met constraints on average and that minimized the empirical average of the objective function, sampled at least 450 times at each point. To optimize performance given the relatively high number of dimensions in our grid searches (Additional file 1: Table S1), we used a grid that dynamically self-refined around the parameter regions in which the objective function was lowest. We used custom-written software that used the Java remote method invocation to distribute jobs to ~1000 single-threaded workers provided by a cluster of 64-core nodes, and dynamically adjusted the grid using aggregated results. Using this setup, a six-dimensional optimization takes ~1–2 days to complete.
As a second step, we used Markov chain Monte Carlo [90] to establish a posterior distribution on the parameters that did at least as well as the best parameter identified by the grid search. Each chain was initiated using that parameter. Burn-in was calculated post hoc so that the autocorrelation of all output parameters decayed to at least 1/e. Each iteration computed an empirical average for the objective function and for constraints using 450 samples; proposed moves were rejected if the empirical average of the objective function was ε higher than the grid-search optimum or if constraints were violated with corresponding functions ε higher than for the grid-search optimum (ε = 0.1 for cell production constraints, ε = 1.0 for fecundity constraints, and ε = 0.01 for the fitness metric unless otherwise specified in Additional file 1: Table S1). Each chain ran for 10,000 iterations. To establish posterior distributions, we thresholded samples along the chain path to keep those that met constraints and did at least as well as the starting point in terms of minimizing the objective function (note that some posterior distributions only contain a large number of repeats of the same point). After thresholding, each optimization had at least 50 samples used to construct the posterior distributions and 95 % credible intervals. We used a parallelized version of our simulations for fast computation of empirical averages at each point, and used custom Python software to drive the process (details of that software will be reported elsewhere); each chain takes ~3 days to complete on a 64-core computer. Detailed results are shown in Additional file 2: Figure S7.
Fitness function for joint optimization of mutation rate and growth rate
We defined a fitness function that captures the effects of delays in the reproductive schedules due to slow germ-cell cycling and of long-term mutational load stemming from replication-dependent mutations. Consider a wild-type population of worms with exponential growth rate r
0 and with gonads that produce progeny with pedigree depth p
0, and a mutant population with growth rate r and gonads that produce progeny with pedigree depth p. Then, making the approximation that all mutations are dependent on replication (see below for a discussion of this assumption), the change in deleterious mutation rate is ΔU = U (p – p
0)/p
0. For selfing species, the selection coefficient for a trait that changes mutation rate by ΔU is s
U
= −ΔU/2 [7, 91, 92]. For an exponentially growing population, the per-generation selection coefficient corresponding to a change in growth rate can be expressed as ln(1 + s
G
) = (r – r
0)/r
0 × ln (N), where N is the number of progeny per generation, and where r = ln(d), where d is the dominant eigenvalue of the population transition matrix resulting from the reproductive schedule and the assumption of a constant speed of embryonic and larval development. Assuming independence of the effects on generation rate and mutation accumulation, the fitness of the mutation population will be 1 + s
U
+ s
G
. We ran optimizations with two numerical values of U. The lower value U = 0.03 reported by [43] was derived from the decrease in fitness of mutation accumulation lines. The higher value U = 0.48 reported by [44] was derived from mutation rates measured by sequencing and the ratio of synonymous to non-synonymous substitution rates. Given our overall conclusion that speed of reproduction takes precedence over pedigree-depth minimization, and given that considering the fraction of mutations that are not dependent on DNA replication would decrease the relative weight of pedigree-depth minimization (although likely not by much, since DNA replication likely plays a preponderant role in mutation accumulation, e.g. [93]), our assumption for the purposes of this computation that all mutations are dependent on DNA replication is conservative.
