Association with pathogenic bacteria affects life-history traits and population growth in Caenorhabditis elegans

Determining the relationship between individual life-history traits and population dynamics is an essential step to understand and predict natural selection. Model organisms that can be conveniently studied experimentally at both levels are invaluable to test the rich body of theoretical literature in this area. The nematode Caenorhabditis elegans, despite being a well-established workhorse in genetics, has only recently received attention from ecologists and evolutionary biologists, especially with respect to its association with pathogenic bacteria. In order to start filling the gap between the two areas, we conducted a series of experiments aiming at measuring life-history traits as well as population growth of C. elegans in response to three different bacterial strains: Escherichia coli OP50, Salmonella enterica Typhimurium, and Pseudomonas aeruginosa PAO1. Whereas previous studies had established that the latter two reduced the survival of nematodes feeding on them compared to E. coli OP50, we report for the first time an enhancement in reproductive success and population growth for worms feeding on S. enterica Typhimurium. Furthermore, we used an age-specific population dynamic model, parameterized using individual life-history assays, to successfully predict the growth of populations over three generations. This study paves the way for more detailed and quantitative experimental investigation of the ecology and evolution of C. elegans and the bacteria it interacts with, which could improve our understanding of the fate of opportunistic pathogens in the environment.


Supplementary Methods
In the following sections, we use t to denote the time of the relevant experiment, and a the age of a worm from the egg-laying time. In the survival and reproduction assays and in the population growth assay, experiments were started (t = 0) when arrested L1 larvae were transferred to bacterial lawns. First, we parametrize a hazard function (death rate) and a reproduction function with respect to the experimental time t. Second, we determine the effective initial age a(t = 0) by comparing the time to adult stage for synchronised larvae and for eggs measured in developmental assays. Third, we define the age-specific birth rate β(a) and death rate µ(a) by adding the corresponding time shift. Let n i be the number of worms observed to be alive at time t i , starting at t 0 . In the absence of censoring, the Kaplan-Meyer estimator of the survival function is equal to the proportion of worms still alive at a given time:

Parametric Survival analysis
Parametric hazard function For the development of the population dynamic model, where the agespecific death rate corresponds to the hazard function, we fit h(t) to the survival data, assuming that it follows a standard mathematical function. In line with previous work (see main text), we chose a Gompertz model, which assumes that the hazard increases exponentially with age: h(t) = β γ e β t , with scale β > 0 and shape γ, which leads to survival S(t) = exp{−γ(e βt − 1)} and to the pdf f (t) = β γ exp{γ + β t − γ e β t }.
Maximum-likelihood fitting of survival functions Given a survival function S(t) with a set of parameter Θ = β, γ, the likelihood of recording a death at time With daily observations, the likelihood of recorded times of death The goodness-of-fit was assessed in two ways: first by running a discrete Cramer-von Mies test using the dgof package (Arnold and Emerson 2011), second by Monte Carlo simulations of observed times of death drawn randomly from the fitted Gompertz distributions.

Developmental assays
In addition to the developmental assay (D egg ) presented in the main text (Fig. 3), which was started from eggs, we ran a similar assay (D L ) starting from synchronised larvae. Arrested L1 larvae were generated as described in the main Methods, and transferred individually to bacterial lawns (t=0). Their status was recorded every two hours from t=0 to t=26 h, and again from t=34 to t=48 h. For each source of food, we then estimated the effective age a 0 of synchronised larvae at t=0 as the difference between the median time to Adult stage t A (as per D L ) and the median age a A at which egg cohorts reached the adult stage (as per D egg ). For example, if the median time from synchronised L1 to Adult was t A = 30 h and the median time age at Adult stage was a A = 36 h, we would infer that synchronised larvae were effectively 6 hour-old at the start of the experiment.

Timing of egg-laying
In the life-history experiments, fecundity at time t was measured as the number of viable larvae from eggs laid between t-1 and t. The fecundity of cohorts was pooled across all adult worms on a plate that were alive at time t-1. As explained in the main text, the recorded number of viable offspring (LRS) was occasionally reduced by early death of mothers. However, the population dynamic model requires that we use the potential LRS that would be achieved in the absence of maternal death. This variable, which we denote as R, was computed for each experimental group by assuming that worms that died between day t and t + 1 achieved on average half of their potential fecundity on that day.
We model egg-laying as a continuous stochastic process, assuming egg-laying events are i.i.d. with probability density function g(t). The predicted proportion of the LRS laid on day t is For a given parametric model g(t; Θ), and daily reproductive success time-series x 1 , ..., x k from a single worm, the likelihood is given by a multinomial distribution, hence the expression for the log-likelihood: The same applies to cohorts of 10 or 25 worms, considering a cohort as a single "super-worm" with a larger number of eggs. Assuming the worms and plates in a given experiment are i.i.d, we can just sum the individual log-likelihoods.
For convenience we assumed that age-specific egg-laying followed a Gamma distribution from the time worms become adults. We used the median value of t A as defined in the previous section, and imposed g(t) = 0∀t < t A . In order to determine fecundity at age a, we use the parametrized function g(t) at time t = a − a 0 , where a 0 is the effective age at t = 0 as per the previous section. More specifically, the age-specific birth rate β(a) was obtained by multiplying g(a − a 0 ) by the potential reproductive success R. We parametrized the model for each experimental group using the average of R (Table 1) to produce Fig. 4 in the main text. Below we also present results from simulations obtained by varying R in line with data from the life-history assays.

Population dynamic model
Let ρ(t, a) be the age-density function of the population at time t: by definition a2 a1 ρ(t, a)da represents the number of individuals in the population with ages in [a 1 , a 2 ] at time t. The total population size at time t is Let β(a) = Rg(a − a 0 ) be the age-specific birth rate, with R and g as defined above. Then with boundary conditions: ρ(t, 0) = B(t), ∀t > 0; ρ(0, a) = φ(a), ∀a > 0, where the initial age distribution φ(a) is given by the experimental design. All experiments were started at time t = 0 by transfering arrested L1 larvae on food. Thus, we can assume a narrow initial age distribution around a L1 .
The birth and death rate functions β(a) and µ(a) were chosen in line with analysis of individual and cohort dynamics. Life-history experiments were initiated at t = 0 with synchronised L1 larvae, hence a L1 < a(t = 0) < a L2 . Ifĥ(t) is the fitted hazard function for survival, we must use death rate µ(a) =ĥ(a − a 0 ), ∀a > a 0 ; the survival of eggs is unknown. The fecundity was measured by counting viable offspring, so it accounts any egg mortality. Therefore we can assume in the model that µ(a) = 0, ∀a < a L1 . In addition, we impose β(a) = 0, ∀a < a adult . The partial differential equation was then solved numerically using the ode.1D function from the deSolve package in R.   Figure B1. Each panel shows the daily proportion of worms alive in each experimental group (black lines), the fitted Gompertz model (red curves), and the 2.5 and 97.5 percentiles of simulated distributions of lifespans from the fitted Gompertz survival functions for each group, using the same number of observations as in the experimental data (pink lines).     Figure B5. Sensitivity to development timing. As in Figure 4, the boxplots represent the data from the population growth experiment, and the green crosses show the predicted numbers of worms in each stage from the model parameterised with data from cohorts of 10 worms. In order to illustrate one of the factors that may have contributed to the relatively poor prediction of the numbers of L4 by the model (particularly in the E. coli and P. aeruginosa groups), we carried out a simple sensitivity analysis to the timing of the onset and duration of the L4 stage. In line with the range of developmental timings reported on Fig. 3, we varied the onset of L4 stage by +/-12h around the mean value, and the average duration was multiplied by a factor ranging from 0.5 to 2. The green vertical bars show the range of the predicted numbers of worms in each stage when both parameters were varied simultaneously.