Microbial phenotypic heterogeneity in response to a metabolic toxin: Continuous, dynamically shifting distribution of formaldehyde tolerance in Methylobacterium extorquens populations

We modeled the dynamics of M. extorquens populations during exposure to formaldehyde with a partial differential equation:

The population is structured by phenotype, with N(x,t) denoting the concentration of cells (CFU*mL^-1) with formaldehyde tolerance x (mM) at time t (hours). The formaldehyde tolerance of a cell is defined as the maximum concentration of formaldehyde in which the cell can grow. The model tracks cells in a well-mixed, closed population as they grow on substrate c at per capita rate r_c (h^-1), die at per capita rate H(x,F), and change phenotype (Fig 9). We assume that phenotypic transitions potentially occur due to two processes, diffusion with coefficient D (mM²*h^-1) and advection with rate v (mM*h^-1) (where mM refers to tolerance level). The death function H(x,F) describes per capita death rate of cells as a function of formaldehyde concentration and is given by:

For each of the three experimental conditions separately, we used maximum likelihood to fit the parameters α, b, D, and v to the data, and used a likelihood ratio test on the nested models to determine the best model structure. For the selection scenario, we began with a 1-parameter model with α > 0 and the other parameters equal to 0, and tested 2-, 3-, and 4-parameter models sequentially, at each step choosing the model with the highest likelihood as long as it was significantly better than the simpler model (results in Table 2). For regrowth (where there is no death), we fit only v and D. In each of the three experimental conditions, our model was able to reproduce the experimental observations extremely well (pseudo-R² = 0.973, 0.993, 0.991 for the formaldehyde selection, methanol regrowth, and succinate regrowth conditions respectively) (Table 2).

Phenotype transition processes change according to growth conditions

We examined not only what the best-fit value was for each of the four parameters of interest, but also whether there was support (by likelihood ratio test) for including each of the parameters. The three experimental conditions differed from one another in both considerations, suggesting that the rate and nature of phenotype transition processes change depending on environment.

In the formaldehyde selection regime, the best model included both death parameters (α = 0.202±0.022 h^-1, b = 0.770±0.147), but for the phenotype transition parameters, only diffusion (D = 0.019±0.006 mM²*h^-1) and not advection. This indicates that the changes we observed in the formaldehyde tolerance distribution during formaldehyde exposure are due not only to death, but also involve phenotype shifts consistent with diffusion. Diffusion leads to the spread of phenotypes consistent with what was observed in the model after 20 hours, including the presence of population density at x = 2 and x = 3 mM (which otherwise would have decreased to below detection in the absence of transitions from higher-tolerance phenotypes) and some cells at x>6 (which were not observed in the initial population, but may have resulted from transitions from lower-tolerance phenotypes) (Fig 10).

For the scenarios involving regrowth of a selected population on formaldehyde-free medium, the parameter estimates were markedly different from those in the selection condition, and from each other. For methanol growth, we found support only for very mild advection toward lower tolerance (v = 0.018±0.006 mM*h^-1). For succinate growth, we found support for both advection (v = 0.285±0.026 mM*h^-1) and diffusion (D = 0.033±0.010 mM²* h^-1); the advection term for succinate was an order of magnitude greater than that for methanol, indicating strong shifts in the direction of lower tolerance when high-tolerance cells are grown on succinate. This agrees with our earlier qualitative observations that growth on succinate, but not on methanol, leads M. extorquens populations to lose formaldehyde tolerance rapidly and in general to undergo diversifying phenotype transitions. Furthermore, it supports our hypothesis that formaldehyde tolerance is associated with methylotrophic growth, and implies that further work is needed not only to understand the mechanism of phenotypic formaldehyde tolerance shifts in M. extorquens but also the regulation of those processes.

The inclusion of b in the model specifically allowed us to test the possibility that a cell's formaldehyde tolerance might determine not only the threshold concentration above which it dies, but also the rate at which it dies (S8 Fig). That is, if b = 0, the death rate is equivalent for all cells regardless of tolerance level as long as tolerance is below the threshold; but for 0<b≤1, cells with lower tolerance levels die more quickly than those with higher tolerance, and the strength of this dependence upon x scales with the value of b. The value of b that best fit our data was 0.770 (Table 2), indicating that formaldehyde-tolerant cells may receive some protection from the tolerance phenotype even at concentrations above their MBC. However, both the experimental results and the model simulation showed a bimodal phenotype distribution after formaldehyde exposure, consistent with a weak dependence of death rate from tolerance level (S8 Fig). Future experiments with a greater number of observations at low-tolerance phenotypes would better elucidate the relationship between tolerance and death dynamics.

Discussion

We have described here a novel example of phenotypic tolerance to a metabolic toxin that is continuously distributed across individual cells in a clonal microbial population, where the phenotypic heterogeneity is present even when all cells inhabit the same environment, but its distribution undergoes dynamic shifts when the population experiences different environmental conditions. Wild-type populations of genetically identical Methylobacterium extorquens cells grown in well-mixed liquid medium contain individuals with maximum formaldehyde tolerance levels ranging from 0 mM to 6 mM (and tolerance at up to 8 mM has been observed after selection); although the distribution is continuous and exponentially decreasing, individuals show binary growth/non-growth phenotypes at any given formaldehyde concentration, and may transition among phenotypic states through both bidirectional phenotype diversification and responsive, directed phenotype shifts.

When might M. extorquens experience formaldehyde in its natural environment, and at what concentrations? M. extorquens excretes formaldehyde during the first stages of the switch between multi-carbon and single-carbon metabolism [27]; in addition, formaldehyde is a metabolic intermediate in the consumption of many lignin-derived aromatic compounds [46,47], and we have observed lignin degraders to excrete formaldehyde into the growth medium at millimolar levels during growth in batch liquid culture on methoxylated aromatic compounds [48]. Thus, it is possible that, in the environment, formaldehyde concentrations in the millimolar range might accumulate transiently on the microscale, especially within cell aggregates such as those observed on plant leaves [49]. However, it is also possible that formaldehyde tolerance is an outward manifestation of a cellular state that has quite a different relevance in the native environment of M. extorquens: high tolerance to local (external) formaldehyde may indicate a high capacity for tolerating (internal) formaldehyde generated through methylotrophic metabolism. As such, this phenomenon may provide insight into processes that are general to many organisms whose central metabolic pathways involve toxic intermediates [50,51].

The mechanisms by which high-tolerance cells can maintain normal growth in the presence of seemingly lethal concentrations of formaldehyde, and by which that ability is transmitted to progeny, remain yet to be elucidated. Our observation of a distinct gene expression profile in the tolerant subpopulation, including the upregulation of stress response genes (universal stress proteins, small heatshock protein 20, and cytochrome bd oxidase) and genes with various roles (transport, regulation, signaling, hypothetical) not known to be linked to formaldehyde, suggest that phenotypic formaldehyde tolerance may involve a system of cellular stress response to damaged, misfolded proteins, and a novel mechanism of formaldehyde tolerance. However, substantial further work will be necessary to uncover the precise mechanism by which tolerant cells are able to maintain normal growth in the presence of formaldehyde, and whether this mechanism is a cause or consequence of the differences in gene expression we observed between the tolerant and sensitive cells. Quantifying gene expression in populations tolerant to a range of different formaldehyde concentrations could shed light onto the relationship between gene expression levels and the continuous nature of the tolerance phenotype distribution.

A further clue to the potential mechanism of tolerance may arise from the association between maintenance of formaldehyde tolerance and methanol growth, and the potentially lower fitness of tolerant cells during succinate growth. These observations hint at a connection between formaldehyde tolerance and methylotrophic metabolism—a connection that would not be evident from the results of gene expression analysis alone. One potential explanation for such a relationship might lie in a physiological feedback loop enabling two different outcomes from initially very small variations in the activity of C₁ pathway enzymes, such as those of the H₄MPT pathway of formaldehyde oxidation, which play a crucial role in removing cellular formaldehyde [30]. Previous reports have found phenotypic heterogeneity in the related strain M. extorquens AM1 in growth rate, cell size, gene expression levels, and ability to switch between carbon substrates [8,52,53]. Assuming that formaldehyde transport into the cytoplasm is diffusion-driven (as no means of active transport has yet been discovered), phenotypic diversity in formaldehyde oxidation capacity would result in a range of internal concentrations across cells. Furthermore, if cells that experience formaldehyde-mediated damage to proteins and other macromolecules begin to lose their capacity to oxidize formaldehyde, this would generate a positive feedback circuit at the protein level that could ultimately determine a binary outcome to whether a cell lives or dies. The fact that we observed no elongation at all from the cells that could not grow at 2.5 mM suggests that this postulated positive feedback mechanism acts on a very fast timescale (i.e., minutes, not hours). This mechanism could explain both the continuous relationship between concentration and abundance of tolerant cells, and the binary nature of the "live normally or die immediately" distribution observed at each specific formaldehyde level, similar to the topology of interactions observed in the case of chloramphenicol resistance described earlier [24].

Similar examples of bi-stable phenotypic outcomes have been described: for instance, pyrimidine-mediated colony type switching in Bacillus subtilis [54], and chloramphenicol resistance [24]. In these cases, the balance between the concentrations of multiple intracellular components determines the phenotype of a cell, with cells pushed away from the threshold into one of two states by a positive feedback loop, and often maintained there by hysteresis (as we observed in the maintenance of formaldehyde tolerance among cells growing on methanol, Fig 8). Making changes that alter the ratios of the components, through control of either gene expression or of the environment, alters the frequency of cells near the threshold and thereby the ratio of cells in each of the two bi-stable states (as we observed when altering formaldehyde concentrations in the medium). Furthermore, these systems are analogous to threshold traits observed in animals, in which observed discontinuous phenotypes (e.g., two phenotypic morphs) arise from a continuously distributed underlying trait [55]. The threshold model of quantitative genetics describes a constant threshold that determines which values of the underlying trait generate each of the phenotypic morphs; evolutionary and environmental changes can lead to shifts in the underlying trait distributions that result in different ratios of the two morphs being observed in different environments. As we have shown here, the continuously distributed trait of M. extorquens formaldehyde tolerance similarly changes with environment; however, in this case the threshold (environmental formaldehyde concentration, which separates tolerant cells from sensitive) is external and may also be manipulated. M. extorquens therefore provides a convenient model system in which to probe the dynamics of threshold traits.

Both stress response proteins and metabolic enzymes are inherited during cell division, so either pathway provides a hypothesis for the heritability of formaldehyde tolerance that we observed after formaldehyde was removed. Lineage dependence has been observed in numerous cases of phenotypic heterogeneity, and is often cited as evidence that the phenotype of interest is dictated by a heritable component of the cell with a moderately long lifetime (e.g., pole age [8], stress-protective protein aggregates [56]). In these cases, the random transitions in phenotype that we modeled as diffusion can be explained by asymmetric partitioning of cell contents [57], whereas the directed transitions that we modeled as advection could be due to responsive up- or down-regulation of the production/degradation or (in)activation of the molecules involved.

The balance between stochastic and responsive phenotype differentiation processes, as mechanisms for increasing population fitness in unpredictable environments, has often been discussed as a question of the relative costs and benefits of each. Both mathematical modeling [58,59] and laboratory evolution [19,60] have demonstrated that either mechanism—or both simultaneously—can be selected for, with evolutionary outcomes dictated by the cost of sensing, the timing of environmental change, the reliability of environmental cues, and the spatial structure of the community. The fact that we observed non-directed general diversification in our populations in some growth conditions, and environment-responsive, directed phenotype shifts in others, strongly suggests that that differentiation in formaldehyde tolerance in M. extorquens is not due simply to unavoidable molecular noise, but rather it is a regulated process conferring a fitness advantage.

If it exists, the nature of such a fitness advantage remains unclear; it is tempting to ask why M. extorquens does not maintain a population composed solely of high-tolerance cells. The evolutionary basis for phenotypic heterogeneity in genetically identical populations of microorganisms is frequently ascribed to diversifying bet-hedging, in which a species in an unpredictably changing environment constitutively generates progeny with multiple phenotypes to ensure that at least a few will thrive in any circumstances [12,61]. The surface of the plant leaf, the native environment of M. extorquens, is indeed unpredictable: cells depend for growth upon a combination of gaseous methanol excreted from plant stomata [62] and other metabolites such as simple organic acids produced by the plant or by other microorganisms [63]. The emission of methanol, dependent on the metabolic state of the plant host and the conductivity of the stomata, undergoes large temporal variation [64]. It is therefore plausible that this environment could select for heterogeneity in formaldehyde toxicity response and in growth on non-methylotrophic substrates. An alternative evolutionary explanation is division of labor, which is an advantageous strategy when a particular activity is beneficial to the population but incurs some cost to the individual carrying it out [59]. In our experiments, the tolerant subpopulation is capable of detoxifying formaldehyde from the growth medium; although in batch culture this occurs too late to prevent the sensitive cells from dying (Fig 2), this might occur differently in a spatially-structured community such as cell aggregates on the surface of a plant leaf.

Importantly, both bet-hedging and division of labor assume fitness tradeoffs among the phenotypes being diversified. Although tolerant cells may be slower than sensitive cells at growth on a multicarbon substrate (S8 Fig), it is unclear whether this disadvantage is substantial enough to explain the low frequency of high-tolerance cells in our M. extorquens populations. Further mathematical modeling could be used to quantify the effects of these growth tradeoffs and to explore conditions that could lead to the evolution of the steady-state distribution that we have observed in the lab. Furthermore, while phenotypic heterogeneity has been studied in great depth both in laboratory populations and in populations simulated by mathematical modeling, there remains a dearth of research on organisms in the natural environments in which they evolved [61,65]. Future experiments examining the dynamics of phenotypic formaldehyde tolerance among M. extorquens cells growing on plant leaves, and among related methylotrophs in other environmental niches, will be essential to establishing the environmental relevance of this phenomenon.

Materials and methods

Bacterial strains and culture conditions

All experiments were conducted with M. extorquens PA1 CM2730, an otherwise wild-type strain that contains a deletion of the cellulose synthesis operon to prevent cell clumping [66]. All cultures were grown at 30 °C in MPIPES mineral medium (30 mM PIPES, 1.45 mM K₂HPO₄, 1.88 mM NaH₂PO₄, 0.5 mM MgCl₂, 5.0 mM (NH₄)₂SO₄, 0.02 mM CaCl₂, 45.3 μM Na₃C₆H₅O₇, 1.2 μM ZnSO₄, 1.02 μM MnCl₂, 17.8 μM FeSO₄, 2 μM (NH₄)₆Mo₇O₂₄, 1 μM CuSO₄, 2 μM CoCl₂, 0.338 μM Na₂WO₄, pH 6.7) [66], either in liquid form or with the addition of 15 g/L Bacto-Agar (BD Diagnostics) for solid medium. Growth substrates consisted either of methanol provided at 15 mM (in liquid medium) or 125 mM (in solid medium); or succinate provided at lower concentrations to ensure the same molarity of carbon as in methanol conditions: 3.5 mM succinate (in liquid) or 15 mM (in solid). Unless otherwise noted, all liquid cultures were grown as a 20 mL volume in 50 mL-capacity conical flasks, capped with Suba-Seal rubber septa (Sigma-Aldrich) to prevent the escape of volatile compounds, with shaking at 250 rpm. Extensive testing in our lab has shown that cultures grown thus are not limited by oxygen availability and grow at the same rate as in flasks with loose lids. Prior to the beginning of each experiment, cultures were streaked from freezer stock onto solid medium and allowed to grow for 4 days to form colonies, and for each biological replicate a single colony was used to inoculate 5 mL of liquid medium in a culture tube and allowed to grow 24 hours. This overnight culture, containing ~2x10⁸ CFU/mL at stationary phase, was used as the inoculum for the growth experiment.

Formaldehyde exposure experiments

Growth and tolerance distribution dynamics in the presence of low concentrations of formaldehyde were assessed by inoculating stationary-phase culture at a 1:64 dilution (for an initial population of ~3x10⁶ CFU/mL) into fresh MPIPES + methanol medium containing formaldehyde at the specified concentration (4 mM unless otherwise stated). Cultures were grown under the conditions described above (20 mL MPIPES, Suba-Seal septa, shaking at 30 °C) and sampled periodically to measure cell viability and formaldehyde concentration. Sampling was conducted as follows: 100 μL of culture was removed through the septum using a sterile syringe and transferred to a trUVue low-volume cuvette (Bio-Rad) to read optical density at 600 nm using a SmartSpec Plus spectrophotometer (Bio-Rad). An additional 200 μL of culture was removed and transferred to a microcentrifuge tube. The culture was centrifuged for 1 minute at 14,000 x g; the supernatant was removed and saved for formaldehyde measurement (see below) and the cell pellet was resuspended in 200 μL fresh MPIPES medium without carbon substrate or formaldehyde. Cells were then subjected to serial 1:10 dilutions in MPIPES to a final dilution of 10⁻⁶. From each of the seven dilutions, three replicates of 10 μL were pipetted onto culture plate containing MPIPES-methanol solid medium to form spots (total: 21 spots per sample per plate type). In experiments measuring formaldehyde tolerance distribution, multiple plate types were used, each type containing a different concentration of formaldehyde (see “Formaldehyde tolerance distributions,” below). The spots were allowed to dry briefly in a laminar flow hood, then lids were replaced and plates were stored in plastic bags and incubated at 30 °C for 4 days before colonies were counted. For each replicate set of seven spots, the two highest-dilution spots with countable colonies were enumerated and summed, then multiplied by 1.1 times the lower of the two dilution factors to calculate the original number of colony-forming units (CFU) in the sample. For each sample, the mean and standard deviation of the three replicate spot series was calculated.

For the measurement of death rates in the presence of high concentrations of formaldehyde, cells were grown in MPIPES medium with methanol and no formaldehyde until they reached an OD₆₀₀ of 0.1 (~ 1x10⁸ CFU/mL, mid-exponential phase), then formaldehyde was added to the desired final concentration and immediately mixed well. Growth conditions were as described above and samples were taken every 20 minutes to measure formaldehyde concentrations and cell viability. Growth rates were measured by fitting a linear relationship between time and the binary logarithm of CFU/mL using the lm function; differences among growth rates were assessed using either the anova function or the t.test function, in the stats package.

Fresh 1 M formaldehyde stock was made weekly by combining 0.3 g paraformaldehyde powder (Sigma Aldrich), 9.95 mL ultrapure water, and 50 μL 10 N NaOH solution in a sealed tube, and immersing in a boiling water bath for 20 minutes to depolymerize. Formaldehyde was measured using the method of Nash [67]. In brief, equal volumes of sample (or standard) and Nash Reagent B (2 M ammonium acetate, 50 mM glacial acetic acid, 20 mM acetylacetone) were combined in a microcentrifuge tube and incubated for 6 minutes at 60 °C. Absorbance was read on a spectrophotometer at 412 nm. For experiments involving large numbers of samples, the same assay was conducted in a 96-well polystyrene flat-bottom culture plate (Olympus Plastics): 100 μL each of sample and Reagent B were combined in each well of a culture plate and the plate was incubated at 60 °C for 10 minutes before reading absorbance at 432 nm on a Wallac 1420 Victor2 multilabel plate reader (Perkin Elmer). A clean plate was used for each assay; each plate contained each sample in triplicate as well as a standard curve run in triplicate.

Genome resequencing of formaldehyde-tolerant subpopulation

Whole-genome sequencing was conducted on an M. extorquens population selected for high formaldehyde tolerance, to determine whether there were any genetic changes associated with the formaldehyde tolerance phenotype. In addition, the genomes of two non-selected M. extorquens populations were prepared, sequenced, and analyzed at the same time using the same methods, as controls to enable us to distinguish mutations specific to formaldehyde tolerance from any that may have accrued in our laboratory strain since the sequencing of the published genome.

Tolerant cells were selected via a 4 mM formaldehyde exposure experiment as described above, and allowed to grow until stationary phase; the tolerance distribution of each population was assayed to confirm that cells were indeed 100% tolerant to 4 mM formaldehyde (see “Formaldehyde tolerance distributions,” below). From each of 3 replicate populations, 2 mL of culture (~2.2x10⁸ CFU/sample) were harvested. Genomic DNA was extracted using the Epicentre Masterpure Complete DNA and RNA Purification Kit (Epicentre/Illumina), following the manufacturer protocol for DNA purification, and the three populations were pooled. Library preparation was carried out by the IBEST Genomic Resources Center at the University of Idaho (Moscow, ID): genomic DNA was used to make shotgun libraries using the TruSeq PCR-free Library Prep Kit (Illumina) with the HiSeq-length insert option (short), amplified using KAPA (Illumina), cleaned using magnetic beads (MagBio), and quantified using fluorometry prior to pooling. The pooled library was quality-checked by Fragment Analyzer (Advanced Analytical) and quantified with the KAPA Library Quantification Kit for ABI Prism (Kapa Biosystems). Sequencing was conducted in a 1x100 run on a HiSeq 4000 at the University of Oregon Genomics and Cell Characterization Core Facility (Eugene, OR); reads were demultiplexed by the facility using bcl2fastq Conversion Software (Illumina).

Genomic data was analyzed for evidence of mutations using breseq 0.32.1 [68] with default settings, using the published genome of M. extorquens PA1 (NC_010172 [69]) as a reference, and comparing with a similar dataset from an M. extorquens CM2730 population that had been grown without formaldehyde. No predicted mutations were found in the formaldehyde-tolerant population that were not also found in the non-formaldehyde population. Three genomic loci (MEXT_RS13110, MEXT_RS12285/ MEXT_RS12290 intergenic region, MEXT_RS02695) carried SNPs identified as "marginal mutation predictions" at a frequency of <33%, but further analysis revealed these to be due to assembly errors in extremely repeat-rich regions that appear in both the tolerant and control samples. Areas containing marginal mutation predictions were further checked by PCR amplification and Sanger sequencing of the original DNA. Primers used were 5'- CTCTCCGCCGAAGTGGT-3' and 5'- GCCTTCCTCGGGTTCAAGGG-3' (for MEXT_RS02695); 5'- CAGGGAACGCTCGTAGAGG-3' and 5'- CCACCGTGAAACGCACCGTA-3' (for MEXT_RS12285-RS12290); and 5'-GTAGACCGCCTCCGAGACTT-3' and 5'- GTAGACCGCCTCCGAGACTT-3' (for MEXT_RS13110). Genome resequencing data have been deposited in the NCBI SRA database under BioProject PRJNA504295.

Cell proliferation assay

Growth versus non-growth phenotypes in bulk liquid culture were assessed using a cell proliferation assay with the non-toxic green fluorescent membrane linker dye PKH67 (Sigma-Aldrich). To ensure that Methylobacterium cells could be easily distinguished from background events in flow cytometry, experiments were conducted using M. extorquens CM3839, a strain identical to CM2730 but constitutively expressing the red fluorescent protein mCherry at the hpt locus [70]. We first conducted a 4 mM formaldehyde exposure experiment using CM3839 to ensure that its population dynamics upon formaldehyde exposure were similar to those of CM2730 (S9 Fig). For the cell proliferation assay, stationary-phase cultures were stained and washed following the manufacturer protocol for PKH67, modified only in that all centrifugation steps were carried out for 1 minute at 14,000 x g. OD₆₀₀ was measured after staining to account for any loss of cells, and inoculation density was adjusted to ensure an initial population of ~3x10⁶ CFU/mL. Stained cells were used in a formaldehyde exposure experiment as described above, with treatments at 0, 1, 2, 3, 4, 5, and 20 mM formaldehyde. Unstained cells were grown alongside the stained cells as a control to assess whether staining affected growth rate or viability; no measurable difference in growth was detected. Samples were taken periodically by syringe and washed of formaldehyde as described above; instead of plating for viability, they were resuspended in 1 M (3%) formaldehyde as a fixative and stored at 4 °C until analysis by flow cytometry. Immediately prior to analysis, cells were centrifuged once more and resuspended in fresh medium to remove excess formaldehyde.

Flow cytometry was conducted in the IBEST Imaging Core at the University of Idaho using a CytoFLEX S benchtop flow cytometer (Beckman Coulter); each sample was analyzed at a flow rate of 10 μL/min for 3 minutes to ensure that an equal volume was examined from each. Output was gated to allow only events with a mCherry signal (ECD-area channel, excitation: 488 nm, emission: 610/620 bp) >10³. Per-cell membrane fluorescence was measured in the FITC-area channel (excitation: 488 nm; emission: 525/40 bp). Data analysis was conducted in R v.3.4.3 [71] using the flowCore package [72] to interpret fcs files and the ggplot2 package [73] to generate plots.

Formaldehyde tolerance distributions

The distribution of formaldehyde tolerance phenotypes within a population was assessed by plating cell cultures onto agar plates containing formaldehyde. MPIPES medium was prepared with agar, autoclaved, and cooled to 50 °C; then methanol and formaldehyde to the desired final concentration were rapidly mixed in, and the agar was poured into 100 mm petri dishes. The dish lids were immediately replaced and plates were cooled on the benchtop. Plates were stored at 4 °C and used within 1 week of pouring. CFU were plated and enumerated as described above. This method has a limit of detection of 1.65x10⁻⁷ (an abundance of 34 CFU/mL is necessary to observe 1 cell per 30 μL plated, and the total cell population tested was 2x10⁸ CFU/mL; therefore the least-abundant subpopulation that could be detected, disregarding the effects of Poisson distributions at lower λ, is one with an average frequency of 1.65x10⁻⁷ within the total population).

To assess whether formaldehyde concentrations in agar culture plates changed over time due to volatilization, plates containing 0, 2, 4, 6, 8, and 10 mM of formaldehyde were assayed for formaldehyde content before and after a 3-day incubation, stored together in the same bag, at 30 °C. A small amount of agar (~0.1 g) was excised from the plate, melted, diluted 1:10 in MPIPES medium, and assayed using the Nash protocol (described above). Each plate was assayed in triplicate. No change in concentration was detected in any of the plates (S10 Fig).

Colony appearance time and growth rates

To measure the effect of formaldehyde damage on colony time and growth rates, colony growth was monitored using Epson Perfection V600 flatbed photo scanners. Formaldehyde exposure experiments were carried out on M. extorquens as described above, but diluted samples were spread-plated rather than spotted onto culture plates: 100 μL of one dilution was spread on each plate. Lids were lined with sterile black felt to reduce condensation and increase contrast, and plates were placed agar-side-down on scanner beds. Each culture contained precisely 30 mL of agar medium, to ensure uniformity of nutrient supply and hydration status across all plates; temperature probes were included between plates on scanner beds to monitor temperatures for consistency. Scanners were placed in a 30 °C incubator and image acquisition was controlled by a computer running Linux Mint, using a cron job for scheduling and a custom bash script employing the utility scanimage to take images once per hour.

Images were processed using a custom Python 3.5.6 script employing scikit-image v.0.12.1 [74] to identify colonies and measure their areas in pixels. Subsequent data analysis was conducted in R: double colonies and non-colony objects were removed from the dataset, colony appearance time was measured as the first timepoint at which colony area measured greater than 100 pixels, and colony growth rate was measured by fitting a linear relationship between time and the binary logarithm of colony area using the lm function. Variability among growth rates was assessed as the median average deviation (MAD) within the colonies from a single timepoint, using the mad function in the stats package. The effect of formaldehyde exposure time on colony appearance time, colony growth rate, and MAD was calculated using simple linear models using the lm function. Further details of both image analysis and colony growth statistics are included in Supporting Information (S2 Fig).

Gene expression analysis

Cultures were prepared by conducting a 4 mM formaldehyde exposure experiment as described above, in biological triplicate. Cultures were harvested at 4 hours of formaldehyde treatment (stressed sensitive cells); 64 hours of formaldehyde treatment or when OD₆₀₀ reached 0.2–0.3 (tolerant growing cells); and 16 hours with no formaldehyde treatment (sensitive unstressed control). Harvesting proceeded as follows: 10⁹ cells were centrifuged in a Beckman tabletop centrifuge (specs) for 5 minutes 4,000xg at 4 °C and supernatant was removed; cells were washed in cold fresh MPIPES medium (no carbon) to remove excess formaldehyde and centrifuged once more to pellet; supernatant was removed, and the pellets were immediately frozen in liquid nitrogen. Harvest volume was 50 mL of exponential-phase culture for the two growing treatments and 150 mL of formaldehyde-stressed cells.

RNA extraction, cDNA generation, library preparation and sequencing were conducted by Genewiz. Briefly, total RNA was extracted via the RNeasy Plus Mini Kit (Qiagen); rRNA was depleted with the Ribo-Zero rRNA removal kit (Illumina); the quality of resulting RNA samples was assessed with an Agilent 2100 BioAnalyzer and Qubit assay. Once generated, the cDNA library was sequenced on an Illumina HiSeq (2x150bp) and raw data in FASTQ format was provided to investigators.

Data analysis was conducted using an analysis pipeline developed by the University of Minnesota Genomics Center and the Research Informatics Solutions (RIS) group at the University of Minnesota Supercomputing Institute. The pipeline uses FastQC [75] to assess the quality of the sequencing data, Trimmomatic [76] to remove the low quality bases and adapter sequences, HISAT2 [77] to align the curated reads to the M. extorquens PA1 genome [GenBank accession: CP000908.1], Cuffquant and Cufnorm from the Cufflinks package [78] to generate FPKM expression values, and featureCounts from the Rsubread R package [79] to generate raw read counts. DESeq2 was used to convert raw count data to normalized counts and all further statistical analysis were carried out using in R [80]. Differentially expressed genes were identified by comparison to the no formaldehyde treatment cells (ie, reference treatment) and had a log₂ fold change > 1 and a False Discovery Rate (FDR) adjusted p-value < 0.001.

Single-cell time-lapse microscopy

To analyze the formaldehyde-dependent heterogeneous response in lag-phase and elongation rates of M. extorquens, we employed single-cell time-lapse microscopy using a phase-contrast inverted microscope (Leica, DMi8) equipped with an automated stage. For image acquisition, we employed a 60x (PH2, NA: 0.7) magnification objective and a sCMOS camera (Hamamatsu, ORCA-Flash 4). Time-lapse imaging was performed under ambient conditions (regulated at 26 °C ± 1 °C) with a 5-minute period using automated routines (Molecular Devices, Methamorph). All reported experiments were performed once; as our setup can accommodate only one condition at a time, the 0 mM and 2.5 mM formaldehyde experiments were each conducted once on different days.

For imaging with single-cell resolution, 2 μL of 10x diluted stationary-phase cells (approximately 6x10⁵ CFU) grown in liquid MPIPES-methanol without formaldehyde were introduced on 1 mm thick agar pads pre-deposited on a glass coverslip. The culture solution was allowed to dry for approximately 10 minutes within a biosafety cabinet and was subsequently covered with a second coverslip. Individual bacteria were therefore immobilized at the agar-coverslip interface, allowing them to grow in two dimensions. Specific locations (200x150 μm²) were monitored in parallel, selected by evaluating the distance between individual cells at the beginning of each experiment so that expanding micro-colonies would not overlap spatially at later time points. On average, 40 microcolonies per location were imaged, with approximately 8–10 individual cells per microcolony in the final timepoint. ImageJ [81] and manual curation were employed for cell segmentation and tracking (S11 Fig).

To prepare the pads, MPIPES-agar medium was prepared as described above ("Bacterial strains and culture conditions"), and melted at 70 °C for approximately 2 hours using a convection oven. Subsequently, methanol was introduced to the agar medium to a final concentration of 125 mM. Specific to the formaldehyde tolerance experiments, formaldehyde was also added to a final concentration of 2.5 mM. In all experiments, a polydimethylsiloxane (PDMS, Sylgard 184, Dow Corning) frame was employed to cover all free edges of the agar pads. This step was employed to minimize formaldehyde evaporation during time-lapse imaging. To fabricate these membranes, PDMS monomer was mixed with its catalyst at a 10:1 ratio, degassed, and cured at 70 °C for approximately 2 hours. Subsequently, the cured PDMS was cut to an area of 25×50 mm², including a 21×47 mm² internal aperture, where the melted agar was subsequently introduced.

Analysis of cell doubling times and microcolony lag times was conducted in R. The wilcox.test function was used to compare the formaldehyde-treated population with the control population by Mann-Whitney Wilcoxon test [82]. The relative effects of formaldehyde treatment and colony lineage were assessed by permutational multivariate ANOVA using distance matrices (PERMANOVA) [83], with the adonis function in the vegan package [84].

For all microscopy experiments, an aliquot of the same culture used for the microscopy experiment was simultaneously grown in batch liquid culture with same concentration of formaldehyde, and CFUs in both the total and formaldehyde-tolerant populations were tracked over time, as a control. In all cases, the growth of the liquid culture proceeded as expected and the frequency of formaldehyde-tolerant cells matched that observed in the microscopy experiment.

Assays on selected formaldehyde-tolerant subpopulations: Fitness trade-offs and formaldehyde-free regrowth experiments

For characterization of the formaldehyde-tolerant subpopulation, a 4 mM formaldehyde exposure experiment was conducted as described above to select for cells with a minimum tolerance level of 4 mM. Once cultures had reached stationary phase (approximately 80 hours), they were used within 4 hours for further experiments.

Assessment of tolerance to antimicrobial drugs and to hydrogen peroxide were carried out using a disk susceptibility test: liquid exponential-phase cultures grown on methanol in standard conditions were diluted to an OD₆₀₀ of 0.3, then spread using a sterile swab onto the surface of MPIPES-methanol-agar plates and allowed to dry. A BBL Sensi-Disc (Becton-Dickinson), impregnated with one antibiotic compound at a set concentration, was then placed on top of the plate. Plates were incubated at 30 °C for 48 hours to allow a lawn to grow; the diameter of the clearing zone around the antibiotic disk was then measured and compared between the naive and selected high-tolerance populations.

To compare growth rates on different carbon substrates, three selected subpopulations, and three naive populations not previously exposed to formaldehyde, were diluted 1:64 into fresh MPIPES medium containing either methanol or succinate, and each of those was transferred into the wells of a 48-well tissue culture plate (Corning Costar) with 640 μL culture in each of three replicate wells. Plates were incubated on a Liconic LPX44 incubator shaker at 650 rpm, and OD₆₀₀ was read using a Wallac 1420 Victor2 plate reader. Growth rates were calculated in R using timepoints during the period of exponential growth (between 5–18 hours for methanol growth and 5–22 hours for succinate), by fitting a linear relationship between time and the binary logarithm of the OD using the lm function in R. Rates were calculated for each well individually; for each biological replicate, the mean and standard deviation of the three replicate wells was calculated. Statistical tests of differences among growth rates consisted of analysis of variance using the anova function from the stats package, and planned contrasts calculated using the lsmeans and contrast functions from the lsmeans package.

To assess tolerance distribution dynamics during formaldehyde-free regrowth, three selected subpopulations were each diluted 1:64 into two batches of fresh MPIPES liquid medium, one batch containing methanol and the other succinate. Dilutions were grown in standard flask conditions as described above, but with loose caps instead of Suba Seals. Cultures were sampled every 4 hours for serial dilution and plating as described above, with each sample plated onto each of 7 plates containing different levels of formaldehyde (0, 2, 4, 6, 8, 10, or 12 mM). Plates were incubated and colonies were counted as described above.

Mathematical model

Growth, death, and phenotype transitions of M. extorquens populations were modeled using a partial differential equation (PDE; Eq 1). The full code for parameter fitting and model selection is included as an R notebook in the Supporting Information, S1 File. Over the time periods for which we ran the model simulation, there was no effective change in the concentration of either of the growth substrates (methanol, succinate; each of which have extremely low half-saturation constants [85,86]) or formaldehyde (Fig 2); these compounds were therefore not explicitly included as time-dependent variables in the model. Furthermore, given that tolerance either has no effect upon growth on methanol, or a mild effect upon succinate growth, we use a single value of r_c for all values of x.

Modeling was conducted in R. The PDE was solved numerically by vectorized ODEs, where each ODE corresponded to a discrete bin within tolerance space of 0.01 mM formaldehyde. A finite difference grid was created using the setup.grid.1D function, and advection and diffusion were calculated using the advection.1D or tran.1D function as appropriate for the different models, in the ReacTran package v. 1.4.3.1 [87]. The vectorized ODEs were solved using the ode.1D function from the package deSolve v1.21 [88], with the lsoda method [89]. Zero-flux boundary conditions for Eq 1 are given by: ∂_xN(0,t) = 0 and ∂_xN(L,t) = 0. The lower boundary was set as x = 0 because tolerance cannot fall below 0. The upper boundary (x = L) for each growth condition was set higher than the highest experimentally observed value, sufficiently large so as not to constrain the upward transition of cells in phenotypic space.

Initial conditions

The experimentally measured distributions of formaldehyde tolerances in the original populations were used as initial conditions in all model runs. Because we found that populations exhibit a slight shift toward higher average tolerance within 2 hours of transitioning from the stationary (non-growing) to the exponential (growing) phase (Fig 7), and because our model specifically focuses on growth phenomena and does not capture behavior in stationary phase, we chose the 2-hour tolerance distribution as our initial condition (S7 Fig) in order to avoid artefacts that would result from using the 0-hour stationary-phase distribution. Three biological replicates were generated and the average of the three used to generate the distribution.

Because the model is continuous, whereas the experimental data were obtained at a resolution of 1.0 or 2.0 mM formaldehyde in tolerance space, the initial conditions for the model were generated by fitting a monotone cubic spline using the function splinefun in R with method “hyman,” to interpolate cell abundances for values of x at intervals of 0.01 mM (S7 Fig). Subsequently, for comparison of model results with experimental results, the model results were summed in 1 mM or 2 mM intervals to re-obtain coarse bins matching the resolution of the experimental data (as in Fig 10).

In the model, x denotes tolerance as the maximum tolerance level of a cell or subpopulation. However, our experimental CFU counts are cumulative in that regard (i.e., the cells that form colonies on culture plates containing 3 mM formaldehyde include those with a maximum tolerance of not just 3 mM but also 4 mM, 5 mM, and 6 mM). For the purposes of carrying out the model, we therefore transformed the empirical tolerance distributions in the experimental datasets to calculate the actual (non-cumulative) number of CFU at each phenotype level using the following formula:

A further adjustment to data from all timepoints was made to account for the experimental limit of detection (34 CFU/mL) in measuring formaldehyde tolerance distributions: growth of originally undetected low-frequency cells in the high-tolerance phenotypes could potentially be mistaken as transition into those phenotypes from lower-tolerance phenotypes. To assess whether such undetected cells would make a difference to the model results, we generated a set of "extended" CFU counts in which we made the extreme correction of adding 1 CFU to the experimental observations (the equivalent of 90.9 CFU/mL for that replicate, or 30.3 CFU/mL if averaged across the three replicates) at high formaldehyde concentrations where the observation had been 0 CFU. This correction was made to the data from all timepoints, according to the following rules:

• If 0 colonies were observed in all technical replicates, add 1 colony to one of the replicates.

• Only do (1) if either:

  • there is an observation of ≥1 colony at a higher concentration than the one being considered, OR

  • it is the first concentration beyond the last observation of ≥1 colony.

We carried out parameter fitting and model selection for the 4 mM formaldehyde selection experiment using both the original dataset and the extended one (as described in the previous paragraph). Although the use of the extended dataset resulted in minor differences in the estimated parameters (Table 2 and S2 Table, S12 Fig), the same form of the model was favored (the 3-parameter model with α, d, and b), and the pseudo-R² (0.973) was slightly better than for the version using the original dataset (0.970). We interpreted this to mean that the inclusion or omission of rare undetected cells at high tolerance levels makes little difference to any biologically relevant conclusions, but that our correction may allow for a better model fit. We therefore proceeded with only the extended dataset for the models describing formaldehyde-free regrowth, and all parameter and statistical values given in the text are those for the extended data.

Parameter estimation

Growth rates (r_c for substrate c, either methanol or succinate, h^-1) and their standard errors were estimated from fitting experimental data of growth of a wild-type, naive population on either succinate or methanol as a primary carbon source, in the absence of formaldehyde (S7 Fig). Linear regression was used to fit the relationship between cell count and time for the exponential portion of three replicate growth curves, using the lm function in R. The growth rate on methanol was 0.195±0.001 h^-1 and the growth rate on succinate was 0.267±0.005 h^-1, where ± denotes 95% confidence interval.

The death rate (α), dependence of death on tolerance (b), diffusion (D), and advection (ν) parameters were estimated using maximum likelihood. Due to the exponential nature of bacterial growth, and because our data sets contain a number of zeroes, we transformed our data using the hyperbolic arcsine function (asinh(x)=ln(x+1+x2)) [60], which is approximately logarithmic but defined at x = 0. The functions lm and logLik in R were used to calculate log-likelihood under a linear model; optimization was carried out in R using optim with the Nelder-Mead method. Standard errors of the parameters were calculated from the Hessian matrix. Fitted values for D and v are given in Table 2 and S2 Table. In addition, we generated a phenotype-independent estimate of α using death rate data from a separate set of time-kill experiments (such as those in Fig 1) at formaldehyde concentrations between 3 and 20 mM. Time-kill curves were conducted as described above, the death rate at each concentration (h^-1) was calculated using linear regression, and the relationship between concentration and death rate was then also calculated using linear regression. By this method, we estimated α as 0.189±0.010 h^-1*mM^-1, which fell within the range of values predicted by fitting α for the models we tested, and within the 95% confidence interval of the estimate of α in the best model. This phenotype-independent estimate was not used in modeling, but helped to verify that the values we obtained by parameter fitting were reasonable.

Model evaluation

For each experimental condition, we used a likelihood ratio test on the nested models using a forward, stepwise procedure to choose the model that best fit the experimental data. For the formaldehyde selection scenario, the "absolute death" model (with α as the only parameter) was used as the null model, and each 2-parameter model (with alpha and either b, D, or v) was compared against it. Of the three 2-parameter models, we chose the one with the highest likelihood as long as it was significantly better than the null model. If a 2-parameter model was chosen, it became the null model and the procedure was repeated to determine whether a 3-parameter model was supported, and subsequently the 4-parameter model. For the regrowth scenarios, where no death due to formaldehyde is possible, we omitted α and b. We first compared two 1-parameter models (v only or D only) against a null model containing neither, and then evaluated whether adding the second parameter was significantly better. LR was calculated as

Supporting information

S1 Fig [pdf]
Formaldehyde concentrations of ≤5 mM allow growth of . at a normal rate, but only after a period of lag; higher concentrations lead to longer lag times.

S2 Fig [pdf]
Image processing pipeline to generate colony growth data from formaldehyde-exposed cultures.

S3 Fig [pdf]
Formaldehyde tolerance may be associated with lower fitness on a multicarbon substrate.

S4 Fig [pdf]
Formaldehyde tolerance distributions in populations are robust across experimental replicates, but vary depending on growth conditions.

S5 Fig [pdf]
Cell proliferation assays support the hypothesis that growth of . in the presence of formaldehyde is due to a small subpopulation of tolerant cells, and that the abundance of tolerant cells decreases with increasing formaldehyde.

S6 Fig [pdf]
The distribution of formaldehyde tolerance within an . population changes over time depending on growth conditions.

S7 Fig [pdf]
Estimation of growth rates and initial conditions for use in the mathematical model.

S8 Fig [pdf]
The parameter (dependence of death rate on formaldehyde tolerance) determines the shape of the population's phenotypic tolerance distribution after exposure to formaldehyde.

S9 Fig [pdf]
Cells expressing mCherry show the same formaldehyde tolerance heterogeneity as wild-type . cells.

S10 Fig [pdf]
Formaldehyde concentrations in agar growth medium are stable over time and reflective of similar concentrations in liquid medium.

S11 Fig [pdf]
Time-lapse microscopy: Cell segmentation and tracking.

S12 Fig [pdf]
Models using extended and original tolerance distributions perform similarly.

S1 Table [pdf]
Tolerant subpopulation shows no difference in sensitivity to antibiotics or hydrogen peroxide.

S2 Table [pdf]
Results of model selection using original data set for fitting (distribution not extended to account for experimental limit of detection).

S1 File [html]
Modeling phenotypic switching in : R notebook.

S1 Data [xlsx]
Excel file with original data shown in each of Figs , , , , and .

S2 Data [csv]
Comma-separated values files with data shown in and used in fitting the model—Required by .

S3 Data [csv]
Comma-separated values files with data shown in and used in fitting the model—Required by .

S4 Data [zip]
Zip file containing full flow cytometry results from experiment shown in .