Weather stressors correlate with Escherichia coli and Salmonella enterica persister formation rates in the phyllosphere: a mathematical modeling study

Case study

The experimental setup for the field studies that provided the bacterial population and weather data used here was previously described by Belias et al. [9]. Briefly, baby spinach and lettuce plants were spray-inoculated with E. coli and S. enterica (Salmonella) onto field plots established in Davis, CA (University of California, Plant Sciences Field Research Facility); Freeville, NY (Homer C. Thompson Research Farm, Cornell University); and Murcia, Spain (La Matanza Research Farm). The spinach and lettuce varieties were selected based on their suitability for baby leaf production: lettuce var. Tamarindo, and spinach var. Acadia F1 and Seaside F1. Four replicate trials at different times of the regional growing season were carried out per location. The plants were spray-inoculated with a 10⁴ CFU/mL cocktail of rifampin-resistant strains of commensal E. coli and attenuated S. enterica serovar Typhimurium (Salmonella), and samples were collected for bacterial cell quantification by plate counts on selective and differential media at 0, 4, 8, 24, 48, 72 and 96 h post-inoculation. Concurrent with leaf sample collection, weather variables (temperature, relative humidity (RH), solar radiation intensity, and wind velocity) were recorded hourly for the respective field locations. The hourly dew point (DP) was calculated as a function of both the hourly temperature and RH.

Model for persister formation on plants

Mathematical modeling to characterize the switch rate from a non-persister bacterial cell (hereafter termed “normal cell”) to a persister cell in the phyllosphere under laboratory conditions was performed as described in our previously published study [24]. Briefly, persister cell fractions were quantified in culturable EcO157 populations after inoculation onto young lettuce plants cultivated in plant growth chambers. Persister cells recovered from the lettuce phyllosphere were identified using the antibiotic lysing method [23]. The greatest persister fraction in the EcO157 population on lettuce in our laboratory investigation above was observed during population decline on leaf surfaces of plants left to dry after inoculation. Using mathematical modeling, we calculated the switch rate from an EcO157 normal to persister cell on dry lettuce plants based on these data [24]. Importantly, our laboratory conditions mimicked inoculation conditions in which E. coli arrived via water on leaves, the surfaces of which progressively dried like under prevailing weather conditions in the field.

Based on the main dynamic observed in the field study data [9] and building on our previous study [24], we assumed that the total enteric pathogen population is composed of (i) non-persister (normal) cells consisting of two sub-populations, characterized by fast (n₁) (CFU/100g) and slow (n₂) (CFU/100g) decay, and (ii) the persister population, leading to the following model from Munther et al. [24]:

where (theta _{n_i})(1/h) is the death rate of the normal cells (subscript i = 1 for fast and i = 2 for slow), (hat p) (CFU/100 g) represents the persister cell population at time t (h), (mu _{hat p}) (1/h) reflects the persister population inactivation rate, α_d (1/h) is the switch rate from normal to persister state, β_d (1/h) is the switch rate from persister to the normal state, and σ ∈ (0,1) is a constant, describing the fraction of persister cells switching back to the normal, slowly decaying state. Equation (1a) and (1b) reflect the assumption that times between switching states are exponentially distributed, using the expected values (frac{1}{{alpha _d}}) (h) and (frac{1}{{beta _d}}) (h) of the respective distributions.

Lacking data for potential persister populations from the field trials, we assumed the persister population is a fraction 1 > k > 0 of the tail population, as observed in Munther et al. [24]. Regarding the model above, this implies that (hat p approx kn_2) for (t ge t^ ast), where (t^ ast approx frac{1}{{theta _{n_1}}}) (the time scale of survival for the fast-decaying population (n₁)). In accord with bi-phasic decay, for (t ge t^ ast), the main dynamics for slow decaying population (n₂) is dictated by (- theta _{n_2}n_2) in Eq. (1b). This suggests that the effective switch rates from n₂ to (hat p) and from (hat p) back to n₂ balance, so that (beta _dsigma hat p approx alpha _dn_2) in Eq. (1b). Following these ideas, we simplified the model in Eq. (1a)–(1d) to:

where we ignored (beta _dleft( {1 – sigma } right)hat p) in (1a) since the decay rate ((theta _{n_1})) dominates. Also, by setting (theta _{hat p} = mu _{hat p} + beta _d(1 – sigma )), and using (beta _dsigma hat p approx alpha _dn_2), we obtained Eq. (2c). Furthermore, because (hat p approx kn_2) for (t ge t^ ast), (theta _{hat p} approx) (theta _{n_2}).

In particular, the assumption that (hat p approx kn_2) for (t ge t^ ast) characterizes the switch rate from normal to persister cells, α_d, as (alpha _d approx kalpha), where α is a hypothetical switch rate assuming that the population is composed only of fast decaying normal cells (n₁) and a hypothetical persister cell population (p). In this case, the hypothetical population p starts small at (widehat {p_0}), initially increases due to switching from population n₁ and then slowly decays as the n₁ population is effectively inactivated (i.e., the tail of the total population is comprised entirely of p). From this perspective we utilized the following equations:

For mathematical justification regarding the relationship (alpha _d approx kalpha), please see the appendix (Supplementary Information).

The utility of the relationship (alpha _d approx kalpha), is twofold. First, we used model fitting (Eqs. (3a)–(3c)) to determine α from the respective field study data [9]. Note that using Eqs. (3a)–(3c), we actually fit for (theta _{n_1}), θ_p, and α using the field study data [9]. Please reference the “model fitting procedure” section as well as the appendix for details concerning the unique determination of the aforementioned parameters, i.e., the practical identifiability of these parameters, and justification regarding the legitimacy of measured tail populations relative to the respective field trial data [9]. Second, because we wanted to examine Spearman’s correlations (corr) between α_d and various weather factors, given a particular weather factor (vec w) across trials (i = 1, ldots ,n), let k be the maximum persister fraction (of the tail) across these n trials, that is, for each i, we have (alpha _{d_i} approx k_ialpha _i), so (alpha _{d_i} lesssim kalpha _i). Thus kα_i represents the maximum persister switch rate for each trial i, and since corr((kvec alpha ,vec w)) =corr((vec alpha ,vec w)), we conducted the correlation analysis with the fitted α values in lieu of the actual persister switch rate α_d.

The assumptions behind our approach are summarized below:

1. A.

   The tails of pathogen populations surviving on plants in the field study [9] are comprised of some fraction k ∈(0,1) of persister cells since their decay rate is quite small and they remain culturable.

2. B.

   Because (alpha _d approx kalpha), we hereafter utilize α from model (3a)–(3c) as the representative persister switch rate.

3. C.

   Given that the experimental context [24] for modeling persister switching occurred during population decline, we only employed trials from Belias et al. [9] that exhibited bi-phasic decay. Namely, we did not include trials in which significant bacterial growth was observed at the time scale of successive data points (the time scale in the field study is on the order of 4–16 h for the 1st day and then 24 h thereafter.)

4. D.

   The switch rate from normal to persister cell is on average a monotonic function of some measure of environmental stress.

Based on assumptions A–D above, we applied the model (3a)–(3c) to published pathogen population size and weather data from four replicate trials in Spain, two in California, and one in NY [9]. More specifically, we fit model (3a)–(3c) to the respective population data in order to:

1. 1.

   determine values for the maximum switch rate α relative to the produce/bacteria type at the field scale,

2. 2.

   describe the correlative relationship between α and weather factors in the respective field trials.

Model fitting procedure

In model (3a)–(3c) above, we supposed dp/dt_t = 0 > 0, i.e., we assumed that bacteria experience stress from the change in conditions from culture growth and inoculum suspension preparation to those on the plant surface and therefore, that persister formation increases in the phyllosphere immediately following inoculation. The report that EcO157 persister formation increases as early as 1 h after inoculation into leaf wash water [23], which could be considered as a proxy for the average oligotrophic environment that bacterial cells experience after spray inoculation onto leaves or through irrigation in the field, supports this assumption. To avoid identifiability issues between the initial persister population (widehat {p_0}) and α regarding the model fits above, we assumed that (widehat {p_0})= 1 ((widehat {p_0}) = 0 gives the same results). Thus, the initial persister population at inoculation is at its lowest, an assumption supported by Munther et al. [24], who observed an average fraction of EcO157 persisters of 0.0043% in the inoculum population. This imparts the largest possible switch rate, α, onto the population, corresponding to the largest and hence most conservative food safety risk.

Let y_k (CFU/100 g of produce) be the average bacteria population measurement at time t_k (h) and let P_k,X (CFU/100 g of produce) represent the model prediction (total population) at time t_k relative to the parameter vector (X = [ {theta _{n_1} , theta_p , alpha } ]^T). Following Eqs. (3a) and (3b), this means that ({{{{{{{mathrm{P}}}}}}}}_{k,X} = n_1left( {t_k,X} right) + p(t_k,X)). Since the population data spans multiple orders of magnitude, we calculated the residuals as (e_{k,X} = log _{10}y_k – log _{10}P_{k,X}). To determine the optimal model fit (see the appendix for details regarding a priori bounds on parameter ranges), we utilized the fminsearch function in MATLAB (MATLAB 2020b, The MathWorks, Inc., Natick, Massachusetts, United States) to determine the parameter vector X that minimizes the 2-norm of the following function F:

Correlation analysis

To provide a statistical foundation from which to relate the switch rate α and measured weather factors, we utilized Spearman and partial Spearman correlation. First, we calculated the Spearman correlation coefficients between α and each of the respective factors: 8-h average of temperature, RH, solar radiation, wind speed post-inoculation, and then we calculated the partial Spearman correlation coefficients for each respective weather factor, while controlling for the other three factors and simultaneously controlling for produce type (using lettuce =1 and spinach =0) (For details regarding why 8-h weather variables were used, see the “model fitting” subsection of the results.) The correlation coefficients were determined using the corr and partialcorr functions in MATLAB 2020b (The MathWorks, Inc., Natick, MA, USA). Considering the significant association of Salmonella α with RH and temperature, we also examined the correlation between α and dew point. Figure 1 presents a logical flow of the statistical analysis. Partial correlations with a P value of less than 0.05 were deemed significant. If the 8-h average of a weather factor exhibited a significant correlation with the switch rate, the 8-h minimum and range of the weather factor were also tested.

Source: Ecology - nature.com