Tunable phenotypic variability through an autoregulatory alternative sigma factor circuit
Genetically identical individuals in bacterial populations can display significant phenotypic variability. This variability can be functional, for example by allowing a fraction of stress prepared cells to survive an otherwise lethal stress. The optimal fraction of stress prepared cells depends on environmental conditions. However, how bacterial populations modulate their level of phenotypic variability remains unclear. Here we show that the alternative sigma factor σV circuit in Bacillus subtilis generates functional phenotypic variability that can be tuned by stress level, environmental history and genetic perturbations. Using single-cell time-lapse microscopy and microfluidics, we find the fraction of cells that immediately activate σV under lysozyme stress depends on stress level and on a transcriptional memory of previous stress. Iteration between model and experiment reveals that this tunability can be explained by the autoregulatory feedback structure of the sigV operon. As predicted by the model, genetic perturbations to the operon also modulate the response variability. The conserved sigma-anti-sigma autoregulation motif is thus a simple mechanism for bacterial populations to modulate their heterogeneity based on their environment.
A combination of single-cell imaging and mathematical modelling reveals a simple mechanism for bacterial populations to modulate their heterogeneity based on their environment and environmental history.
Cells live in a changeable environment and experience a wide range of environmental stresses. Bacterial populations have evolved strategies to survive these stresses. One strategy is for all cells to immediately respond to stress with the activation of the relevant stress response pathway (Hilker et al, 2016). Alternatively, a bacterial population can generate a broad range of cellular states, which allows it to hedge its bets against the changeable environment (Veening et al, 2008b). Noise in gene expression has been proposed as a mechanism for generating phenotypic variability in genetically identical cells (Raj & van Oudenaarden, 2008; Martins & Locke, 2015). This phenotypic variability has also been shown to be affected by changes in the cellular environment, such as a shift in stress level or growth conditions (Megerle et al, 2008; de Jong et al, 2012; Mitosch et al, 2019), as well as by previous ‘priming’ stresses (Mitosch et al, 2017). However, how the bacterial population regulates individual cell decisions to modulate the fraction of cells that enter an alternative transcriptional state remains unclear (Fig 1A).
The σV mediated lysozyme stress response pathway in Bacillus subtilis is an ideal model system to examine how bacterial populations can tune their phenotypic variability. σV is an extracytoplasmic function (ECF) alternative sigma factor. Alternative sigma factors replace the ‘housekeeping’ sigma factor, σA, in the RNA polymerase holoenzyme and redirect it to regulons that control distinct regulatory programmes. They have already been shown to display a high level of gene expression variability in B. subtilis (Locke et al, 2011; Young et al, 2013; Cabeen et al, 2017; Park et al, 2018), and the σV activation pathway is both well characterized and specific to one stress condition, which greatly simplifies analysis of its activation.
σV is the only pathway activated in response to C-type lysozyme (Guariglia-Oropeza & Helmann, 2011; Ho et al, 2011; Ho & Ellermeier, 2012). Lysozyme is produced by animals as part of their innate immune system and kills bacteria by cleaving the peptidoglycan of the cell wall between the N-acetylmuramic acid residue and the β-(1,4)-linked N-acetylglucosamine (Lal & Caplan, 2011). In its inactive form, σV is bound to its anti-sigma factor RsiV, a transmembrane protein (Fig 1B). If lysozyme is present, RsiV binds to lysozyme (Hastie et al, 2014; Hastie et al, 2016) and activates a signal transduction cascade to release σV. First RsiV undergoes a conformational change that allows signalling peptidases to cleave RsiV at site-1 (Hastie et al, 2014; Castro et al, 2018; Lewerke et al, 2018) (Fig 1B). Bacillus subtilis has five type 1 signal peptidases, of which the two major peptidases are SipS and SipT (Tjalsma et al, 1998). Either SipS or SipT is sufficient for site-1 cleavage (Castro et al, 2018; Ho & Ellermeier, 2019). The truncated RsiV can then be cleaved by RasP (a site-2 protease), which results in the release of σV (Hastie et al, 2013; Hastie et al, 2014) (Fig 1B).
Once σV is released it can bind to RNA polymerase and redirect transcription to the σV regulon. The regulon contains genes that allow adaptation to lysozyme stress, such as the penicillin-binding protein PbpX and the d-alanyl-d-alanine carrier protein ligase DltA. Both these genes are also part of other sigma factor pathways (Guariglia-Oropeza & Helmann, 2011; Ho et al, 2011). Like many alternative sigma factors, one of σV's targets is its own operon. The sigV operon includes sigV (the gene that codes for σV) and its anti-sigma factor rsiV, and so its activation results in ‘mixed’ positive and negative feedback loops. In addition, the operon contains oatA, one of the main lysozyme resistance genes and yrhK, which is of unknown function (Hastie & Ellermeier, 2016). OatA contributes to lysozyme adaptation by transferring an acetyl group to the C-6-OH position of N-acetylmuramic acid in the peptidoglycan, thus blocking cell wall cleavage by lysozyme (Bernard et al, 2011; Ho et al, 2011).
Although the molecular mechanisms underlying σV activation by lysozyme have been elucidated, previous work was carried out using bulk experiments that average out single-cell dynamics (Guariglia-Oropeza & Helmann, 2011; Ho et al, 2011) and mask cell-to-cell heterogeneity. This heterogeneity can be crucial in identifying and distinguishing between regulatory strategies (Munsky et al, 2012). In this study, we used single-cell time-lapse microscopy of fluorescent reporters for σV activity to characterize σV induction dynamics in individual cells in response to lysozyme stress. We found that upon induction by sub-lethal levels of lysozyme, σV is activated heterogeneously, with some cells activating σV rapidly, whereas some cell lineages do not activate σV for multiple generations. This heterogeneity is functional, as cells that respond to an initial sub-lethal stress are more likely to survive a subsequent lethal stress of lysozyme. Through experiment and modelling, we found that these dynamics can be explained solely by the autoregulatory feedback of σV on its own operon. Our model predicted that the observed heterogeneity could be tuned by environmental history and by genetic perturbations, which we confirmed experimentally. The conserved sigma-anti-sigma autoregulation motif is thus a simple mechanism for bacterial populations to tune their levels of phenotypic variability.
Individual cells activate σV heterogeneously under lysozyme stress
To characterize σV activation dynamics, we first constructed a B. subtilis reporter strain containing a chromosomally integrated fluorescent reporter for σV activity, PsigV-YFP. This strain also contained a reporter for the housekeeping sigma factor σA (PtrpE-RFP). PtrpE-RFP expression was used as a constitutive control and to aid image analysis (Materials and Methods). We then used time-lapse microscopy to examine σV activity in individual cells grown in the mother machine microfluidic device (Wang et al, 2010). This device allows long-term tracking of hundreds of individual cells trapped at the ends of channels. It also allows for fast switching between media conditions during imaging (Materials and Methods). We first measured σV expression dynamics in response to a step change to a sub-lethal concentration of 1 µg/ml lysozyme (Fig 1C and D, Appendix Fig S1 and Fig EV1 and Movie EV1). The first cells begin to respond to the lysozyme stress with raised PsigV-YFP levels within two frames (20 min), suggesting rapid activation given the ˜ 15 min maturation time of the YFP reporter protein (Appendix Fig S2). However, our time-lapse movies revealed that PsigV-YFP was heterogeneously activated at the single-cell level (Fig 1C). While 20% of cells reached the half maximum of σV activity within 70 min of being exposed to lysozyme, it took 200 min (approximately four generations) for 90% of all cells to activate σV (Fig 1D and Appendix Fig S3B).
To test whether the observed heterogeneous activation of σV in response to lysozyme was due to the growth conditions in the mother machine, we investigated σV activation dynamics in liquid culture (Appendix Fig S4) and using an alternative microfluidic device (Appendix Fig S5). Both conditions showed similar heterogeneous σV activation to that observed in the mother machine. Therefore, the observed heterogeneity in PsigV-YFP is independent of the experimental setup. Previous work has shown that heterogeneous gene expression could be due to intrinsic noise in the expression of the PsigV-YFP reporter (Elowitz et al, 2002), rather than due to heterogeneous σV activity. To test whether the observed heterogeneity was due to the intrinsic variability of the PsigV-YFP reporter, we constructed a strain containing chromosomally integrated PsigV-YFP and PsigV-mTurq reporters and exposed it to 1 µg/ml lysozyme in liquid culture. In these snapshot experiments, the expression of PsigV-YFP and PsigV-mTurq was highly correlated (R2 = 0.85) (Appendix Fig S6). Thus, the observed heterogeneity in PsigV-YFP reflects changes in σv activity and not intrinsic variability of the PsigV-YFP promoter.
Next, we asked whether the observed heterogeneity in PsigV-YFP is modulated by the level of lysozyme applied. We examined PsigV-YFP expression after the application of 0.5, 1, 2 and 4 µg/ml lysozyme. These values were all below the previously reported minimal growth inhibitory concentration of 6.25 µg/ml (Ho et al, 2011) and led to a transient decrease in growth rate (Appendix Fig S7). To measure the distribution of σV activation times, for each time point we calculated the fraction of cells that had crossed the half-maximum of their respective final σV value (meaning the cell had activated σV). We found that when increasing the lysozyme concentration from 0.5 to 4 µg/ml the heterogeneity in σV activity was reduced (Figs 1F and EV1 and Appendix Fig S8). The time for 90% of cells to activate their σV pathway decreased from 300 min (approximately six generations) for 0.5 µg/ml to 100 min (approximately two generations) for 4 µg/ml lysozyme (Appendix Fig S3). At the same time, the fold change in induction between the unstressed σV activity and the steady-state σV activity under lysozyme stress increased from ˜ 190 for 0.5 μg/ml to ˜ 520 for 4 μg/ml (Fig 1G and Appendix Fig S1). We also observed that under a 4 µg/ml concentration of lysozyme some cells (8 and 21% in two different repeat experiments) appeared sick and were wider than usual cells. These cells also overshot their new σV activity steady-state before relaxing to it. We removed these cells from our analysis (Appendix Fig S9), although including them did not affect our results (Appendix Fig S10). We also observed that the fraction of cells activating σV increased with increasing lysozyme in the alternative microfluidic device, although in this device movies were stopped before all cells activated due to crowding of the cells in the chip (Appendix Fig S5F). Taken together, our results reveal that the level of phenotypic variability in σV activation times is tuned by stress levels, with the heterogeneity in σV activation reduced as lysozyme levels increase.
We examined whether the differences in individual cell states before applying a lysozyme stress could predict how a cell responds to the stress. There was no correlation between either growth rate or cell size before lysozyme application and response time after lysozyme application (Appendix Figs S11 and S12). Alternative sigma factors in B. subtilis display gene expression variability even in the absence of stress (Locke et al, 2011; Park et al, 2018). We found that, in the absence of stress, the coefficient of variation of PsigV-YFP was over four times higher than that of the constitutive control (0.65 ± 0.29 vs 0.14 ± 0.05) (Appendix Fig S13). Cells which had elevated levels of PsigV-YFP before the addition of stress were more likely to activate σV instantaneously (Fig EV2 and Appendix Fig S14). However, some cells with low PsigV-YFP levels also activate σV instantaneously (Fig EV2 and Appendix Fig S14). Therefore, while noise in PsigV-YFP expression before stress contributes to the observed heterogeneity in σV activation after the addition of lysozyme, it is not the only cause.
Given the heterogeneity in σV activation times, we examined whether activating σV early had any effect on the survival against future lethal concentrations of lysozyme. Cells that were exposed to 20 μg/ml of lysozyme for 20 min all died within 1 h (Fig 2 and Appendix Fig S15). However, if the cells were first exposed to a priming stress of 1 μg/ml of lysozyme for 30 min, which heterogeneously induced σV, and then subsequently to 20 μg/ml of lysozyme for 20 min, some cells survived the high lysozyme stress (Fig 2A). We chose this priming stress level and duration as previous experiments had shown (Appendix Fig S4) that it ensured heterogeneous activation of σV, with a large fraction of cells not having turned on before the second higher stress. We found that the short exposure to sub-lethal concentrations of lysozyme (during the priming stress) improved survival to subsequent lethal stress levels from 0 to 12.5 ± 2.7% (See Materials and Methods, Fig 2B and Appendix Fig S15). Cells that survived until the end of the movie, 280 min after the 20 μg/ml of lysozyme, had on average 1.57 ± 0.33 fold higher PsigV-YFP levels than perishing cells immediately before the application of the lethal concentration of lysozyme (Fig 2C).
The heterogeneous activation of PsigV-YFP is solely due to σV and its anti-sigma factor RsiV
We next attempted to understand how the single-cell activation dynamics of the σV pathway are generated. First, to test whether the heterogeneity that we observed in σV activation times was due to lysozyme stress activating different stress response pathways, we analysed the genome-wide transcriptional response of cells to 1 μg/ml lysozyme. We carried out RNA-seq 30 min after the addition of stress in the wild type and in the ΔsigV knockout. We chose this stress level and duration as we had seen in previous experiments that it resulted in heterogeneous σV activation (Appendix Fig S4). As previously reported, only the sigV operon was strongly (> 5 fold induction) induced by lysozyme (Guariglia-Oropeza & Helmann, 2011) in the wild type (Fig 3B). The lysozyme resistance gene pbpX, which is part of the sigV regulon, was also upregulated in the WT (4.9 fold induction), but was not upregulated in the ΔsigV background (0.88 fold induction), consistent with the known role of sigV in its activation (Guariglia-Oropeza & Helmann, 2011). The ΔsigV strain did not show any genes with strong (> 5 fold induction) induction (Fig 3C). DltA, which is known to contribute to lysozyme resistance and is part of the sigV regulon (Guariglia-Oropeza & Helmann, 2011), was not upregulated in the WT or the ΔsigV strain. Our results suggest that no other pathway is strongly induced by 1 µg/ml lysozyme. We therefore hypothesized that the observed PsigV-YFP heterogeneity is due to the σV circuit itself (Ho et al, 2011; Hastie et al, 2013; Hastie et al, 2014).
To test which components of the σV circuit play a role in the heterogeneous activation of σV, we constructed IPTG-inducible overexpression constructs of all genes in the sigV operon (sigV, rsiV, oatA, yrhK), as well as genes in the σV circuit (sipS, rasP) and the sigV regulon (pbpX) in a PsigV-YFP reporter background (Fig 3A). We then investigated σV activation dynamics under lysozyme stress after full induction of each circuit component. To do this, we took single-cell snapshots of PsigV-YFP expression of cells grown in liquid culture. The cells were either unstressed or were stressed for 30 min with 1 μg/ml lysozyme. As a metric for the variability in activation time of σV in response to lysozyme, we calculated the fraction of cells with PsigV-YFP levels larger than a threshold above the unstressed level of each strain (See Materials and Methods). Single-cell snapshots of PsigV-YFP expression revealed that overexpressing either yrhK, sipS, rasP or pbpX only had a minor effect on the fraction of activated σV cells (Fig 3D) or on the level of PsigV-YFP expression (Appendix Fig S16). For cells overexpressing sigV itself, no further activation of σV occurred on the addition of lysozyme (Fig 3D). This was because PsigV-YFP expression was already at a high level prior to stress (Appendix Fig S16). Conversely, overexpression of rsiV caused cells not to express PsigV-YFP at all, preventing activation of all cells on addition of lysozyme (Fig 3D and Appendix Fig S16). Finally, overexpressing oatA, which blocks lysozyme cleavage of the peptidoglycan, shuts off the activation of σV (Fig 3D and Appendix Fig S16). However, when we increased the lysozyme concentration to 20 μg/ml in the presence of oatA overexpression, the heterogeneous expression of PsigV-YFP reappeared (Fig EV3), suggesting oatA is not responsible for the heterogeneous activation dynamics. We repeated the experiment for rsiV overexpression, but increased RsiV did not increase protect against lysozyme as a concentration of 20 μg/ml lysozyme killed all cells.
To further validate the importance of sigV and rsiV as compared to oatA for the observed heterogeneity in σV activation, we constructed deletion mutants of sigV, rsiV and oatA. Only the deletion of oatA left the activation σV dynamics unchanged (Fig 3E). In the sigV mutant, PsigV-YFP levels did not increase in response to lysozyme stress (Fig 3C and Appendix Fig S17). Deleting rsiV caused all cells to have high PsigV-YFP expression even before the addition of lysozyme and the addition of lysozyme did not activate the system any further (Appendix Fig S17). These findings suggest that the heterogeneity in σV activation only depends on σV and its anti-sigma factor RsiV (Fig 3F). We found that the leakiness of the inducible sigV construct (Phyperspank-sigV, without addition of IPTG) increased the fraction of activated cells on the addition of lysozyme, as well as causing an increase in the steady-state levels of PsigV-YFP before and after stress (Appendix Figs S18 and S19). The leakiness of the RsiV construct (Phyperspank-rsiV) caused the opposite effect. These effects were still apparent for inducible constructs with lower leakiness (Pspank-sigV, Pspank-rsiV (Appendix Figs S18 and S19)), confirming that sigV activation dynamics are sensitive to these two system components.
Mathematical modelling reveals that PsigV-YFP expression dynamics can be explained by the ‘mixed’ σV autoregulatory feedback loop
Based on our transcriptome and overexpression analysis, we constructed a model of a simplified sigV operon consisting of sigV and its anti-sigma factor rsiV, with positive autoregulation of the operon by active σV (Fig 3F). We did not include oatA, as it is not required for the heterogeneous activation of σV (Fig EV3). RsiV binds σV to form an inactive complex σV-RsiV. On stress application, RsiV is degraded and σV is free to activate the operon. The model was simulated using Gillespie-type stochastic simulations, tracking the changes in copy number of each species of the system by simulating the individual reaction events (Gillespie, 1977). With the exception of the production of σV and RsiV (which were implemented through a Hill function of σV activity with a Hill coefficient of 2), all reactions were modelled using mass action kinetics (See Materials and Methods). A Hill coefficient > 1 was required to generate heterogeneous activation dynamics as a degree of ultrasensitivity is needed in the system to amplify the response to molecule fluctuations. We note that other networks, including a single positive feedback loop, can generate heterogeneous activation dynamics (See Materials and Methods and Appendix Fig S20), but that this is the simplest model that can simulate the role of both σV and RsiV in the regulation of sigV activation.
We searched a range of biologically feasible parameters and were able to find parameters that capture the heterogeneous σV activation in response to lysozyme stress (Figs 4A and EV4). The sigV operon components were assumed to be stable, so the dilution rate was set to approximately match the division rate observed in experiments. Simulations yielded plausible copy numbers for the number of sigma factor molecules (Jishage & Ishihama, 1995; Jishage et al, 1996). In addition, we verified that the heterogeneous activation behaviour is robust to perturbations to the system parameters (Appendix Fig S21). Using our model, we were able to recreate several aspects of the experimental data, including the dependency of induction time (Fig 4B) and steady-state levels of sigV expression on levels of lysozyme stress (Fig 4C). Finally, the heterogeneous activation dynamics were also modulated by small increases in the baseline production rate of either σV or RsiV (Appendix Fig S22), qualitatively matching the effects of the leakiness of the inducible σV or RsiV construct observed in experiment (Appendix Fig S18).
Our model consists of a mixed positive and negative feedback loop. We tested the requirements of this feedback for the dynamics by modelling a feedback-broken system, with constitutive expression of sigV and rsiV. For a range of constitutive expression, the dynamic range of PsigV-YFP expression for the feedback-broken system on addition of lysozyme was less than that of the WT system (Fig EV5A and Appendix Fig S23). This reflected the role of the feedback loop in amplifying the system dynamics. To test this prediction experimentally, we constructed a strain with no autoregulation of the sigV operon by knocking out the sigV operon and replacing it with a sigV operon driven by an inducible promoter. This allowed us to study the system at different steady-state expression levels (by varying IPTG induction level) to the WT system. We found that the fold change induction of the WT on addition of lysozyme was at least 4.5 times higher than that observed in the inducible operon strain, regardless of the IPTG induction level (Fig EV5B and Appendix Fig S24), matching the behaviour observed in our model (Fig EV5A).
Phenotypic variability is tuned by doubling copy numbers of sigV operon genes
To further test the assumptions of our model, we predicted the effect on σV activation dynamics of introducing a second copy of each component of the sigV operon (sigV, rsiV or both sigV and rsiV), with each component driven by the sigV promoter (Fig 5A). Our model predicted that these perturbations would modulate the heterogeneity of the system's response to lysozyme stress. A second copy of rsiV meant that large fluctuations in sigV expression were required to kick the system into the high σV expression state (Fig 5B and Appendix Fig S25). This increased the cell-to-cell variability in response times to lysozyme stress, as compared to the WT in the simulations (Appendix Fig S25). Conversely, with a second copy of sigV, or a second copy of both sigV and rsiV, smaller fluctuations in sigV expression are required to kick the system into the high σV expression state. Cell-to-cell variability in response times to lysozyme therefore decreases in these simulations (Fig 5 and Appendix Fig S25). In fact, a second copy of sigV caused an increase in the levels of σV expression even before the addition of lysozyme in the simulations (Appendix Fig S25).
To test these predictions, we constructed strains containing a second copy of either sigV, rsiV or sigVrsiV driven by the sigV promoter. As predicted by our model, snapshots of sigV expression after induction by 1 μg/ml of lysozyme revealed that a second copy of sigV or sigVrsiV reduced the observed heterogeneity in σV activation, while a second copy of rsiV caused an increase in the heterogeneity in σV activation (Fig 5C and Appendix Fig S26). Finally, adding a second copy of sigV alone caused the activation of σV even in the absence of stress, as predicted. These results validate our model assumptions and also show that the heterogeneity in sigV activation is easily tuned by simple genetic perturbations.
The heterogeneous activation dynamics of σV are dependent on environmental history
Given that the σV activation dynamics appear sensitive to the baseline levels of σV, RsiV and the σV-RsiV complex, any perturbations to these levels should affect σV activation times. We hypothesized that cells should have elevated levels of sigV operon components after a lysozyme stress is removed due to the recent high expression of the operon components. The elevated levels of σV stored in stabilized σV-RsiV complexes should cause all cells to respond immediately to a subsequent reapplication of lysozyme. We investigated this hypothesis in our model. We simulated the application of lysozyme stress, which was then removed for the equivalent duration of several cell cycles. The same level of stress was then reapplied and the heterogeneity in σV activation disappeared (Fig 6A and Appendix Fig S27). However, as the delay before the second application of lysozyme was increased (allowing system components to relax to pre-stress levels), the heterogeneity gradually returned (Fig 6B and C and Appendix Fig S27). Thus, our simulations predict the existence of a temporary molecular memory of the environmental history.
To test whether the σV circuit exhibits such a memory, we grew cells in the mother machine device under 1 μg/ml lysozyme stress before removing the stress for 6 h (approximately seven cell cycles). All cells stopped activating σV as soon as lysozyme was removed, as indicated by the decay in YFP fluorescence (Fig 6D and Movie EV2). The decay in fluorescence was tightly synchronized across the population and YFP decayed with a half-life time of 51 ± 1 min, which was similar to the cell cycle time of 51 ± 13 min. When reapplying a 1 µg/ml lysozyme stress after 6 h (approximately seven generations), all cells responded instantaneously. As predicted by the model, the heterogeneity in σV activation was lost (Fig 6D, Appendix Fig S28A and Movie EV2). Also as predicted by the model, increasing the duration of the break in stress to 12 h (approximately 14 generations) resulted in heterogeneous activation dynamics similar to those observed in response to the first application of lysozyme stress (Fig 6E and F, Appendix Fig S28B and C and Movie EV3). This return to a heterogeneous response reflects the system's loss of transcriptional memory through reduction in the sigV operon component concentrations to pre-stress exposure levels. These results held for two different versions of the mother machine microfluidic device (Appendix Fig S29). Taken together, our results show how the autoregulatory σV circuit can generate heterogeneous activation dynamics that can be tuned by stress level, genetic perturbations and the environmental history of the cell.
Here, we report a general mechanism for a bacterial population to tune its phenotypic variability based on stress levels, genetic architecture and environmental history. Using quantitative single-cell microscopy and microfluidics, we found that the activation of the alternative sigma factor σV in response to lysozyme stress was heterogeneous. While some cells activated their σV pathway immediately, others could take up to six generations to activate their σV pathway. The observed phenotypic variability plays a functional role, as cells that respond to a sub-lethal stress were more likely to survive a subsequent higher stress application (Fig 2). Through experiments and modelling, we found that this heterogeneity could be understood by the ‘mixed’ positive and negative feedback of σV activating both itself and its anti-sigma factor RsiV. Alternative sigma factors are a common regulatory system in prokaryotes, often controlling stress response and virulence pathways. The ‘mixed’ feedback loop is also a common motif, suggesting that this motif can be a general mechanism for bacterial populations to tune their phenotypic variability.
Our modelling and experiments found that recent exposure to lysozyme stress modulates the heterogeneity observed on a second stress application, even though the system turns off after removal of the first lysozyme stress. The key to this behaviour appears to be the ‘mixed’ feedback loop, as it allows amplified levels of inactive σV-RsiV complexes in each cell after stress. These complexes can be cleaved by a subsequent addition of lysozyme, releasing σV to activate its operon. Similar transcriptional memories of previous stress have been observed in bacterial systems, although typically not to modulate phenotypic heterogeneity. For example, other pathways such as the heat stress response in B. subtilis (Runde et al, 2014) or the oxidative stress response in yeast (Kelley & Ideker, 2009) have been shown to have a transcriptional memory of past conditions. Often this transcriptional memory is facilitated through an autoregulatory positive feedback loop that can lock the cell into an ON state that is heritable through cell divisions (Novick & Weiner, 1957; Biggar & Crabtree, 2001; Xiong & Ferrell, 2003; Acar et al, 2005; Hashimoto et al, 2013; Lambert & Kussell, 2014). However, it is difficult for a single positive feedback loop to allow the system to be OFF but primed for future stress, as we find to be the case for the ‘mixed’ feedback loop in the σV pathway. Dilution during growth causes heterogeneous activation of sigV to eventually return when levels of σV -RsiV drop to pre-stress levels, so the memory is qualitatively different from that generated by a positive feedback loop locking a system ON. However, we find that the sigV transcriptional memory remains for several generations. In the future, it will be important to investigate whether the ‘mixed’ feedback loop also tunes the levels of phenotypic diversity by environmental history in other systems. Interestingly, computational studies have proposed that a ‘mixed’ feedback loop structure in the MAR operon in Escherichia coli allows the tuning of the fraction of cells prepared to survive antibiotic exposure (Garcia-Bernardo & Dunlop, 2013). Additionally, the mixed feedback loop mechanism could be compared to other mechanisms proposed to allow the modulation of phenotypic variability, such as multi-site phosphorylation (Libby et al, 2019) or threshold-based mechanisms in toxin–antitoxin modules (Rotem et al, 2010).
While simple, our model allows qualitative matches to data. In future, it will be important to increase the complexity of the model to make more precise predictions of the behaviour of the sigV system. One aspect of the system that can be modelled in more detail is how noise in gene expression is generated in the circuit. In our model, we do not model transcription and translation separately and assume uncorrelated noise for each reaction channel. A more detailed model could involve characterizing the noise in terms of bursts of transcription and translation (Friedman et al, 2006). In turn, this would require experiments to characterize the noise in transcription, such as single-molecule FISH (Raj & Oudenaarden, 2008). Our assumption of uncorrelated noise is also a simplification as, for example, we have modelled the degradation events as uncorrelated, which may not hold as these are primarily caused by dilution. Additionally, the system requires ultrasensitivity (in the form of a Hill coefficient n > 1 in the operon production term) in order to amplify molecule fluctuations. While other sigma factor response systems have been shown to utilize ultrasensitivity (Narula et al, 2012; Narula et al, 2016), there is no known source of it in the σV circuit (neither the binding of σV to RsiV, nor its operon, is cooperative). Further research should examine possible sources of ultrasensitivity in the circuit, one possibility being sigma factor competition for RNA polymerase (Park et al, 2018).
Cells that activate sigV quickly after a priming stress have a higher chance of surviving a subsequent high stress (Fig 2B). This points to a potential benefit of early activation of sigV. Future work should examine the costs of early sigV activation to see if the heterogeneous activation dynamics we observe represent a bet-hedging strategy (Veening et al, 2008a; Veening et al, 2008b), where early responding cells suffer a fitness penalty in return for protection against future stress. We do not observe a growth rate difference in cells that activate sigV earlier compared to later, suggesting that early responders are not suffering a growth rate penalty. However, it is possible that we are missing small growth rate effects, as the time resolution of our mother machine experiments only allows approximately 5 time points to be measured per cell cycle. Our experiments do indicate that high constitutive expression of sigV or oatA reduces the growth rate in bulk culture (Appendix Fig S30). It could also be that the fitness penalty is due to early sigV activation blocking cells from responding to stress with other alternative sigma factors, as it appears alternative sigma factors compete for RNA polymerase (Nyström, 2004; Park et al, 2018). Evolution experiments under changeable stressful environments could reveal whether the heterogeneity in activation and transcriptional memory that we observe evolve to optimally match the external environment.
We have found that a combination of noise in gene expression and a mixed feedback loop can generate tunable phenotypic diversity in an alternative sigma factor circuit. Our work will be of utility for synthetic biology, where alternative sigma factors are a promising system for engineering orthogonal gene circuits (Rhodius et al, 2013; Bervoets et al, 2018; Pinto et al, 2018). Going forward, it will be important to observe whether alternative sigma factors more generally are used as a mechanism for bacterial populations to tune phenotypic diversity. This is particularly the case given that alternative sigma factors often control pathways critical to pathogenicity and resistance to antibiotics. Noise in sigma factor activity has already been observed in multiple alternative sigma factor circuits in B. subtilis (Locke et al, 2011; Narula et al, 2016; Park et al, 2018), as well as for the general stress response sigma factor in E. coli (Patange et al, 2018). Additionally, an alternative sigma factor plays a role in generating phenotypic diversity in Mycobacteria (Sureka et al, 2008; Ghosh et al, 2011). Two other obvious candidates for further study are the pathogens Enterococcus faecalis and Clostridioides difficile (Ho & Ellermeier, 2019), which also have a σV pathway that is responsive to lysozyme.
Materials and Methods
Strains and media
All strains are derivatives of the PY79 background strain (Appendix Table S1). Deletions were generated by replacing genes of interest with an antibiotic resistance cassette by recombination of a linear DNA fragment homologous to the region of interest. All strains had a ΔytvA::neo deletion insertion. YtvA is a blue light sensor, and it was deleted to avoid any activation of the general stress response pathway σB by the microscope illumination (Gaidenko et al, 2006; Locke et al, 2011). For strains used in mother machine experiments, cells were made immotile by inserting a ΔhaG:erm deletion cassette in order to improve the loading into the microfluidic device. Additionally, all strains contained a house keeping σA promoter-driven mCherry for segmentation and as a constitutive control (Locke et al, 2011).
Cells were routinely grown in Spizizen's Minimal Media (SMM) (Spizizen, 1958). It contained 50 µg/ml tryptophan as an amino acid source and 0.5% glucose as a carbon source. Cultures were started from frozen stock in SMM and grown overnight at 30°C to an OD between 0.3 and 0.8. The overnight cultures were resuspended to an OD of 0.01 and regrown to an OD of 0.1 at 37°C.
A Nikon inverted Ti-E microscope (Nikon, Amsterdam, Netherlands) with a Nikon Perfect Focus System (PFS) hardware auto-focus was used. All images were acquired using either a Photometrics CoolSnap HQ2 CCD camera (Photometrics, Tucson, AZ, USA) or a Photometrics Prime sCMOS camera (Photometrics, Tucson, AZ, USA). The microscope stage was enclosed in an incubator (Solent Scientific, Segensworth, UK) which was set to 37°C for all experiments. Illumination was provided by a LED lamp (CoolLED, Andover, UK). Epifluorescence was provided by a Lumencore Solar II light engine (Lumencore, Beaverton, OR, USA). Chroma filters (Chroma, Bellows Falls, USA) #41027 for the RFP channel, #49003 for the YFP channel and #49001 for the CFP channel were used. All experiments were done with a phase 100x Plan Apo (NA 1.4) objective (Nikon, Amsterdam, Netherlands). Metamorph (Molecular Device, Sunnyvale, CA, USA) controlled the camera, the motorized stage (Nikon, Amsterdam, Netherlands) and the microscope.
Agarose pads were made by pipetting 1 ml of 1.5% (wt/vol) low melt agarose (Merck, Darmstadt, Germany) dissolved in PBS onto a 22 mm2 cover glass slide. Immediately after pipetting, another cover glass slide was placed on top of the agarose creating a “sandwich” with the agarose in the middle. Once the agarose had solidified, the top cover glass was removed and the agarose was cut into squares of approximately 5 mm × 5 mm with a scalpel.
Once the cultures had reached an OD of 0.1, they were split into aliquots of 1 ml to which different concentrations of lysozyme from hen eggs white (Sigma Aldrich, St. Louis, MO, USA) were added. The aliquots were then incubated at 37°C for 30 min. For experiments with IPTG-inducible strains, 1 mM of IPTG was added to the aliquots for 60 min before addition of lysozyme, in order to allow for full induction of the promoter before stress. Different concentrations of lysozyme were then applied for 30 min before snapshot measurements. For snapshots, 2 µl of cell culture was pipetted onto the pads, which were left to dry and then laid face down onto a cover glass-bottom dish (#HBSt-5040, WillCO-dish, Amsterdam, Netherlands). The glass dish was sealed with parafilm and put under the microscope to image. Single-cell data were extracted using custom MATLAB (Mathworks, Natick, USA) scripts based on the Schnitzcells package (Young et al, 2011).
Overnight cultures were grown and then resuspended to an OD of 0.01 as described above. Once the cultures had reached an OD of 0.1, they were resuspended to an OD of 0.001 for loading into CellAsic B04 microfluidic chips (Merck, Darmstadt, Germany). Cells were loaded into the microfluidic chips with a pressure of 4–6.5 psi for 2–4 s. The following lysozyme concentrations were investigated with the CellAsic setup: 0, 0.5, 1 and 4 µg/ml. Fresh media was perfused into the chip with a pressure of 1 psi. After 60 min of growth in standard SMM, the media was switched to media containing lysozyme. Cells were imaged at regular intervals (every 10 min), and the acquired movies were analysed with the standard Schnitzcells package for MATLAB (Young et al, 2011).
Sylgard 184 polydimethylsiloxane (PDMS) from Dow Corning (Midland, MI, USA) was used to replica-mould the microfluidic devices. A 10:1 (base:curing agent) PDMS mixture was cast onto the epoxy master and cured overnight at 65°C. In the next morning, the chips were cut out and the inlets were punched using a 0.75 mm biopsy puncher (#504529, WPI, Sarasota, FL, USA). To remove any uncured PDMS from the chips, these were chemically treated in a pentane (Sigma Aldrich, St. Louis, Mo, USA) bath for 90 min and then washed twice for 90 min in acetone (Sigma Aldrich, St. Louis, Mo, USA). The chips were then left to dry overnight in the fume hood and were ready to use the next day (Gruenberger et al, 2013).
The chips were plasma bonded to a glass-bottom dish (HBSt-5040, Willco Wells, Amsterdam Netherlands) by treating the chip and the glass dish surfaces with a plasma (Femto Plasma System, Diener, Ebhausen, Germany) for 12 s at a power of 30 W and a pressure of 0.35 mbar. The chips were then put into the oven for 10 min at 65°C, in order to strengthen the bonding. In the meantime, a bovine serum albumin (BSA) passivation buffer was prepared at a concentration of 20 mg/ml in SMM. The chips were then passivated with this buffer for 1 h at 37°C.
Cells were grown overnight as for the other movie experiments. 10 ml of the regrown cultures with an OD of 0.1 was spun down for 10 min at 3,000 g (4,000 rpm) and resuspended in 0.1 ml to a final OD of 10. This cell suspension was loaded into the growth channels by spinning the chip for 5 min at 3,000 rpm with a spin coater (Polos SPIN150i, SPS, Harlem, the Netherland).
Setup of microfluidic mother machine experiments
Perfusion was provided by a Fluigent pressure pump setup. Briefly, reservoirs were pressurized using a computer controllable pressure controller (MFCS-EZ, 1 bar, Villejuif, Fluigent, France). A rotary selection valve (M-Switch, Fluigent, Villejuif, France) allowed for switching between different reservoirs. The flow rate was set to 1 ml/h and by using a flow sensor (M-Flow Sensor, Fluigent, Villejuif, France) which feeds back onto the pressure controller.
For the experiments on the memory in the σV response, syringe pumps (Legato 200, KD Scientific, Holliston, USA) were used. Switching was provided by an electronic switching valve (MXP9960, Rheodyne, Lake Forest, USA). The flow rate remained at 1 ml/h.
For all shown memory experiments, the JLB221 strain (ΔespH) was used instead of JLB130, unless stated otherwise. This was done in order to prevent clogging of the microfluidic chip during these long-term experiments due to the production of extracellular matrix. The two strains had the same PsigV-YFP dynamics in response to lysozyme (Appendix Figs S31 and S32).
Cells were inoculated in SMM and grown overnight to an OD of 0.3–0.8. In the morning, the cells were resuspended to an OD of 0.01 in SMM + 1 μg/ml lysozyme and grown to an OD of 0.1. Cells were then loaded as described in the “Experimental preparation” section.
In the mother machine, the cells were grown for 4 h with 1 μg/ml lysozyme before removing the stress. The stress was then reapplied after 6 or 12 h.
Growth curve determination
Cultures were started from frozen stock in SMM and grown overnight at 30°C. The overnight cultures were resuspended to an OD of 0.02 and grown for 2 h at 37°C before adding 1 mM of IPTG. The cells were then grown for another hour before adjusting the OD to 0.1 and aliquoting the culture into 2 ml samples (with one 2 ml sample used for each time point), and then, the samples were returned to grow at 37°C. A sample was taken every hour, and the OD600 was determined.
Cumulative activation time
The activation time in mother machine experiments was calculated as the time between switching from SMM to lysozyme and the time point when then PsigV-YFP passed its half maximum (Appendix Fig S1). The cumulative fraction was then calculated as the cumulative fraction of cells with PsigV-YFP values higher than the half maximum of their final values (representing cells that have activated) (Appendix Figs S1 and S8).
Removal of overshooting cells
For 4 μg/ml lysozyme, the PsigV-YFP activity overshot its steady-state activity. Overshooting cells were sick and also wider. Thus, overshooting cells could be removed based on their width. First, the single-cell width traces were smoothed with a Gaussian filter, and their maximum value was determined. All single-cell traces with a maximum width larger than the mean width + 6 sigma of cells before the stress were removed (Appendix Figs S9 and S10).
Calculation of surviving cells in priming experiment
The lysis curve (Fig 2B) was calculated as the fraction of surviving cells 280 min after the addition of 20 µg/ml lysozyme. We only examined the survival rate of the top three cells in each channel at the time of the addition of 20 µg/ml lysozyme to avoid problems with cells being washed out of the chip during the 280 min after the stress. Only one surviving lineage (the longest) was counted from each of these top three cells, to avoid cell division artificially increasing the number of survivors. The survival rate at the end of the movie using this method was 12.5 ± 2.7% (Fig 2B). To verify this method, we also used an alternative approach. In this second method, the end of each channel was approximated using the average cell position before a cell left the channel in the ˜ 100 frames before the stress was added. Cells that at their last recorded position after the addition of 20 µg/ml lysozyme stress were within six standard deviations of the end of the channel were removed from the analysis as they were assumed to have left the channel rather than died. This second approach gave qualitatively similar results, with a survival rate of 13.6 ± 3.1%.
Statistical analysis of the effects of higher PsigV-YFP levels before lysozyme stress application on response time
As the distributions in Appendix Fig S14 were not normal, we use the non-parametric Kolmogorov–Smirnov test (ks test) to verify the null hypothesis that cells with higher PsigV-YFP values before lysozyme stress application would have the same activation time distribution as cells with lower PsigV-YFP expression. Cells with high PsigV-YFP were defined as cells with a YFP fluorescence higher than the mean YFP fluorescence before the addition of lysozyme plus one standard deviation. All other cells were defined as low PsigV-YFP. A P-value below 0.05 was used to reject the null hypothesis. The ks test was performed using the built in MATLAB function kstest2.
We found that for 0.5, 1 and 2 μg/ml lysozyme the null hypothesis was rejected. Thus, cells with higher PsigV-YFP values before lysozyme stress application had a shorter activation time. For 4 μg/ml lysozyme, the null hypothesis was accepted. Therefore, cells with higher PsigV-YFP values before lysozyme stress application had the same activation time as cells with lower PsigV-YFP values.
Calculation of fraction of activated cells from snapshots
Cells which, after the application of lysozyme stress, had a higher mean PsigV-YFP expression than the mean PsigV-YFP expression plus six standard deviations before stress were defined as having activated. The fraction of activated cells was then calculated as the number of cells which had activated normalized by the total number of cells (Figs 3 and 5).
Lysozyme treatment before RNA extraction
Overnight cultures were prepared as for snapshot or movie experiments. Briefly, cultures were started from frozen stock in SMM and grown overnight at 30°C to an OD between 0.3 and 0.8. The overnight cultures were resuspended to an OD of 0.01 and regrown to an OD of 0.1 at 37°C. Once the cultures had reached an OD of 0.1, they were split into aliquots of 10 ml to which either 0 or 1 µg/ml lysozyme from hen eggs white (Sigma Aldrich, St. Louis, MO, USA) were added. The aliquots were then incubated at 37°C for 30 min.
RNA extraction and library preparation
10 ml of B. subtilis at OD 0.1–0.2 was centrifuged in 15-ml falcon tubes for 10 min at 3,000 g (4,000 rpm). Cell pellets were resuspended in 1 ml of Qiagen RNAprotect Bacteria Reagent and flash frozen in liquid nitrogen.
Defrosted cells were pelleted by a quick centrifugation. After removing the supernatant, cells were resuspended in 1 ml buffer RLT from Qiagen RNeasy kit. Resuspended cells were transferred in a Fastprep Lysing Matrix B tube (MP Bio) and processed in Fastprep apparatus 45 s at speed 6.5 M/s. 700 µl of the supernatant, containing lysed cells, was transferred to a new microcentrifuge tube, to which 500 µl of absolute ethanol was added. After vortexing, the lysate was transferred to a RNeasy spin column and centrifuged 15 s at > 9,400 g (10,000 rpm). RNA purification was then carried out following the instructions of the Qiagen RNeasy kit. RNA quality and integrity were assessed on the Agilent 2200 TapeStation, and RNA concentration was assessed using Qubit RNA HS assay kit. Library preparation was performed using ScriptSeq™ Complete Kit (Illumina, BB1224), for 2 µg of high integrity total RNA (RIN > 8). The libraries were sequenced on a NextSeq500 using paired-end sequencing of 75 bp in length.
The raw reads were analysed using a combination of publicly available software and in-house scripts. We first assessed the quality of reads using FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were aligned to the B. subtilis PY79 transcriptome (NCIB no.CP006881) using Salmon (Patro et al, 2017). Read counts for each gene were imported using the tximport R package (Soneson et al, 2015). Genes differentially expressed (P-value < 5%) between the WT and ΔsigV mutant, or in response to lysozyme treatment were identified using the DESeq2 R package (Love et al, 2014). The RNA-seq data of our study can be found here (Gene Expression Omnibus: GSE171761, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171761).
We constructed a mathematical model of the σV circuit that was based on the core components that we found experimentally to modulate the heterogeneity in σV activation (Fig 3). The model consists of σV, RsiV and the σV-RsiV complex. Free σV that is not bound in a complex activates the production of both σV and RsiV. Lysozyme stress was modelled as activating the cleavage of the σV-RsiV complex, releasing σV while degrading RsiV. We note that although heterogeneous activation can be simulated through just a positive feedback loop, (Tiwari et al, 2011; Frigola et al, 2012) (Appendix Fig S20), our model is the simplest model that could be compared to the mutations including both sigV and rsiV (Fig 5).
We used the following set of biochemical reactions to model the dynamics of σV, RsiV and σV-RsiV.
The production rate of σV and RsiV was assumed to follow a Hill function, where the production rate of both σV and RsiV were identical and activated by σV. The production rate is as follows: v0 + v*(σV)ⁿ/((σV)ⁿ + Kⁿ), where v0 (= 0.1 molecules/min, is the operon leakage), v (= 2.5 molecules/min, is the maximal operon activity), K (= 60 molecules, is the apparent dissociation constant) and n (= 2, is the Hill coefficient) are parameters.
We assumed a constant and identical dilution rate for all three components of the system. The dilution rate is set by the parameter kdeg (= 0.01 min−1).
We assumed that σV and RsiV would bind to each other to form a complex and that this complex could dissociate. We also assumed the binding rate to be higher than the dissociation rate. The binding rate is set by the parameter kB (= 10 molecules−1 min−1) and the dissociation rate by the parameter kD (=5 min−1).
We assumed that the σV/RsiV complex would be cleaved by lysozyme, with σV getting released and RsiV degraded. The cleavage rate is set by the product of the two parameters L (which attains variable values, is the amount of lysozyme stress experienced by the system) and kC (= 0.05 min−1, is the base cleavage rate of the complex). Splitting the rate into two parameters allows setting more natural lysozyme input values (such as 1 and 2).
The model was implemented in the Julia programming language using the Catalyst.jl modelling package. Simulations were made using the DifferentialEquations.jl package (Rackauckas & Nie, 2017b). To account for the stochastic nature of the system, we use a Gillespie-type model (Gillespie, 1977). Here, we tracked the copy numbers of the three components of the system (σV, RsiV and σV-RsiV), and their change due to the individual reaction events. We used Gillespie's direct stochastic simulation algorithm to determine the time to, and which, the next reaction event in the simulation should be. We used DifferentialEquations.jl's SSAStepper method to simulate the model. For more details, please see the implementations in the files provided. The model parameters were hand-picked by carrying out a coarse-grained search of feasible values. The degradation rate was set to be similar to what would be produced by the bacteria division rate in the experiments. This yields plausible copy numbers for the number of sigma factor molecules (Jishage & Ishihama, 1995; Jishage et al, 1996). See Appendix Table S2 for parameter values. The heterogeneous activation behaviour is robust with respect to perturbations to the selected parameters (Appendix Fig S21). Finally, we also note that a continuous SDE model, based on the chemical Langevin equations, could reproduce the features of the system, (Appendix Fig S33), suggesting that that behaviour does not depend on choice of model approach (Gillespie, 2000; Rackauckas & Nie, 2017a).
The RNA-seq data produced in this study can be found here:
Gene Expression Omnibus GSE171761: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171761.
Model simulation, data analysis and figure plotting code, as well as source data for main and extended view figures can be found here: https://gitlab.com/slcu/teamJL/schwall_etal_msb_2021.
We thank John Sauls and Professor Suckjoon Jun for the kind gift of a mother machine epoxy mould. We thank Jin Park for help with the mother machine analysis code, Michael Elowitz for sharing plasmids, Craig Ellermeier for sharing strains, critical reading of the manuscript and helpful discussions, Hugo Tavares for help with RNA-seq analysis, Katie Abley for critical reading of the manuscript and Niklas Korsbo for helpful discussions. This research was made possible by the award of a European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement 338060. The work in the Locke laboratory is further supported by a fellowship from the Gatsby Foundation (GAT3272/GLC). Torkel Loman has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 721456.
CPS, TEL and JCWL conceived and designed the study. CPS, TEL, BMCM and JCWL analysed and interpreted the data and wrote the article. CPS, VK and TL performed the experiments. SC did the RNA-seq experiments. CV and TS constructed strains. TEL and BMCM developed the mathematical model.
Conflict of interest
The authors declare that they have no conflict of interest.
- Acar M, Becskei A, Oudenaarden A (2005) Enhancement of cellular memory by reducing stochastic transitions. Nature 435: 228–232CrossrefCASPubMedWeb of Science®Google Scholar
- Bernard E, Rolain T, Courtin P, Guillot A, Langella P, Hols P, Chapot-Chartier M-P (2011) Characterization of O-acetylation of N-acetylglucosamine: a novel structural variation of bacterial peptidoglycan. J Biol Chem 286: 23950–23958CrossrefCASPubMedWeb of Science®Google Scholar
- Bervoets I, Van Brempt M, Van Nerom K, Van Hove B, Maertens J, De Mey M, Charlier D (2018) A sigma factor toolbox for orthogonal gene expression in Escherichia coli. Nucleic Acids Res 46: 2133–2144CrossrefCASPubMedWeb of Science®Google Scholar
- Biggar SR, Crabtree GR (2001) Cell signaling can direct either binary or graded transcriptional responses. EMBO J 20: 3167–3176Wiley Online LibraryCASPubMedWeb of Science®Google Scholar
- Cabeen MT, Russell JR, Paulsson J, Losick R (2017) Use of a microfluidic platform to uncover basic features of energy and environmental stress responses in individual cells of Bacillus subtilis. PLoS Genet 13: e1006901CrossrefPubMedWeb of Science®Google Scholar
- Castro AN, Lewerke LT, Hastie JL, Ellermeier CD (2018) Signal peptidase is necessary and sufficient for site 1 cleavage of RsiV in Bacillus subtilis in response to lysozyme. J Bacteriol 200: e00663-17CrossrefPubMedWeb of Science®Google Scholar
- Eldar A, Chary VK, Xenopoulos P, Fontes ME, Losón OC, Dworkin J, Piggot PJ, Elowitz MB (2009) Partial penetrance facilitates developmental evolution in bacteria. Nature 460: 510–514CrossrefCASPubMedWeb of Science®Google Scholar
- Elowitz MB, Levine AJ, Siggia ED, Swain PS (2002) Stochastic gene expression in a single cell. Science 297: 1183–1186CrossrefCASPubMedWeb of Science®Google Scholar
- Friedman N, Cai L, Xie XS (2006) Linking stochastic dynamics to population distribution: an analytical framework of gene expression. Phys Rev Lett 97: 168302CrossrefCASPubMedWeb of Science®Google Scholar
- Frigola D, Casanellas L, Sancho JM, Ibañes M (2012) Asymmetric stochastic switching driven by intrinsic molecular noise. PLoS One 7: e31407CrossrefCASPubMedWeb of Science®Google Scholar
- Gaidenko TA, Kim T-J, Weigel AL, Brody MS, Price CW (2006) The blue-light receptor YtvA acts in the environmental stress signaling pathway of Bacillus subtilis. J Bacteriol 188: 6387–6395CrossrefCASPubMedWeb of Science®Google Scholar
- Garcia-Bernardo J, Dunlop MJ (2013) Tunable stochastic pulsing in the Escherichia coli multiple antibiotic resistance network from interlinked positive and negative feedback loops. PLoS Comput Biol 9: e1003229CrossrefCASPubMedWeb of Science®Google Scholar
- Ghosh S, Sureka K, Ghosh B, Bose I, Basu J, Kundu M (2011) Phenotypic heterogeneity in mycobacterial stringent response. BMC Syst Biol 5: 18CrossrefCASPubMedWeb of Science®Google Scholar
- Gillespie DT (1977) Exact stochastic simulation of coupled chemical reactions. J Phys Chem 81: 2340–2361CrossrefCASWeb of Science®Google Scholar
- Gillespie DT (2000) The chemical Langevin equation. J Chem Phys 113: 297CrossrefCASWeb of Science®Google Scholar
- Gruenberger A, Probst C, Heyer A, Wiechert W, Frunzke J, Kohlheyer D (2013) Microfluidic picoliter bioreactor for microbial single-cell analysis: fabrication, system setup, and operation. J Vis Exp 50560PubMedWeb of Science®Google Scholar
- Guariglia-Oropeza V, Helmann JD (2011) Bacillus subtilis σV confers lysozyme resistance by activation of two cell wall modification pathways, peptidoglycan O-acetylation and D-alanylation of teichoic acids. J Bacteriol 193: 6223–6232CrossrefCASPubMedWeb of Science®Google Scholar
- Hashimoto M, Seki T, Matsuoka S, Hara H, Asai K, Sadaie Y, Matsumoto K (2013) Induction of extracytoplasmic function sigma factors in Bacillus subtilis cells with defects in lipoteichoic acid synthesis. Microbiology (Reading, Engl) 159: 23–35CrossrefCASPubMedWeb of Science®Google Scholar
- Hastie JL, Williams KB, Ellermeier CD (2013) The activity of σV, an extracytoplasmic function σ factor of Bacillus subtilis, is controlled by regulated proteolysis of the anti-σ factor RsiV. J Bacteriol 195: 3135–3144CrossrefCASPubMedWeb of Science®Google Scholar
- Hastie JL, Williams KB, Sepúlveda C, Houtman JC, Forest KT, Ellermeier CD (2014) Evidence of a bacterial receptor for lysozyme: binding of lysozyme to the anti-σ factor RsiV controls activation of the ecf σ factor σV. PLoS Genet 10: e1004643CrossrefPubMedWeb of Science®Google Scholar
- Hastie JL, Ellermeier CD (2016) Proteolytic activation of extra cytoplasmic function (ECF) σ factors, In Stress and Environmental Regulation of Gene Expression and Adaptation in bacteria FJ Bruijn (ed.), pp 344–351. Hoboken, NJ: John Wiley, Sons Inc.Wiley Online LibraryGoogle Scholar
- Hastie JL, Williams KB, Bohr LL, Houtman JC, Gakhar L, Ellermeier CD (2016) The Anti-sigma factor RsiV Is a bacterial receptor for lysozyme: co-crystal structure determination and demonstration that binding of lysozyme to RsiV is required for σV activation. PLoS Genet 12: e1006287CrossrefPubMedWeb of Science®Google Scholar
- Hilker M, Schwachtje J, Baier M, Balazadeh S, Bäurle I, Geiselhardt S, Hincha DK, Kunze R, Mueller-Roeber B, Rillig MC et al (2016) Priming and memory of stress responses in organisms lacking a nervous system. Biol Rev Camb Philos Soc 91: 1118–1133Wiley Online LibraryPubMedWeb of Science®Google Scholar
- Ho TD, Hastie JL, Intile PJ, Ellermeier CD (2011) The Bacillus subtilis extracytoplasmic function σ factor σV is induced by lysozyme and provides resistance to lysozyme. J Bacteriol 193: 6215–6222CrossrefCASPubMedWeb of Science®Google Scholar
- Ho TD, Ellermeier CD (2012) Extra cytoplasmic function σ factor activation. Curr Opin Microbiol 15: 182–188CrossrefCASPubMedWeb of Science®Google Scholar
- Ho TD, Ellermeier CD (2019) Activation of the extracytoplasmic function σ factor σV by lysozyme. Mol Microbiol 112: 410–419Wiley Online LibraryCASPubMedWeb of Science®Google Scholar
- Jishage M, Ishihama A (1995) Regulation of RNA polymerase sigma subunit synthesis in Escherichia coli: intracellular levels of sigma 70 and sigma 38. J Bacteriol 177: 6832–6835CrossrefCASPubMedWeb of Science®Google Scholar
- Jishage M, Iwata A, Ueda S, Ishihama A (1996) Regulation of RNA polymerase sigma subunit synthesis in Escherichia coli: intracellular levels of four species of sigma subunit under various growth conditions. J Bacteriol 178: 5447–5451CrossrefCASPubMedWeb of Science®Google Scholar
- Jong IG, Veening J-W, Kuipers OP (2012) Single cell analysis of gene expression patterns during carbon starvation in Bacillus subtilis reveals large phenotypic variation. Environ Microbiol 14: 3110–3121Wiley Online LibraryCASPubMedWeb of Science®Google Scholar
- Kelley R, Ideker T (2009) Genome-wide fitness and expression profiling implicate Mga2 in adaptation to hydrogen peroxide. PLoS Genet 5: e1000488CrossrefCASPubMedWeb of Science®Google Scholar
- Lal M, Caplan M (2011) Regulated intramembrane proteolysis: signaling pathways and biological functions. Physiology (Bethesda) 26: 34–44CrossrefCASPubMedWeb of Science®Google Scholar
- Lambert G, Kussell E (2014) Memory and fitness optimization of bacteria under fluctuating environments. PLoS Genet 10: e1004556CrossrefPubMedWeb of Science®Google Scholar
- Lewerke LT, Kies PJ, Müh U, Ellermeier CD (2018) Bacterial sensing: a putative amphipathic helix in RsiV is the switch for activating σV in response to lysozyme. PLoS Genet 14: e1007527CrossrefPubMedWeb of Science®Google Scholar
- Libby EA, Reuveni S, Dworkin J (2019) Multisite phosphorylation drives phenotypic variation in (p)ppGpp synthetase-dependent antibiotic tolerance. Nat Commun 10: 5133CrossrefPubMedWeb of Science®Google Scholar
- Locke JCW, Young JW, Fontes M, Hernández Jiménez MJ, Elowitz MB (2011) Stochastic pulse regulation in bacterial stress response. Science 334: 366–369CrossrefCASPubMedWeb of Science®Google Scholar
- Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15: 550CrossrefCASPubMedWeb of Science®Google Scholar
- Martins BMC, Locke JCW (2015) Microbial individuality: how single-cell heterogeneity enables population level strategies. Curr Opin Microbiol 24: 104–112CrossrefCASPubMedWeb of Science®Google Scholar
- Martins BMC, Tooke AK, Thomas P, Locke JCW (2018) Cell size control driven by the circadian clock and environment in cyanobacteria. Proc Natl Acad Sci USA 115: E11415–E11424CrossrefCASPubMedWeb of Science®Google Scholar
- Megerle JA, Fritz G, Gerland U, Jung K, Rädler JO (2008) Timing and dynamics of single cell gene expression in the arabinose utilization system. Biophys J 95: 2103–2115CrossrefCASPubMedWeb of Science®Google Scholar
- Mitosch K, Rieckh G, Bollenbach T (2017) Noisy response to antibiotic stress predicts subsequent single-cell survival in an acidic environment. Cell Syst. 4: 393–403.e5CrossrefCASPubMedWeb of Science®Google Scholar
- Mitosch K, Rieckh G, Bollenbach T (2019) Temporal order and precision of complex stress responses in individual bacteria. Mol Syst Biol 15: e8470Wiley Online LibraryPubMedWeb of Science®Google Scholar
- Munsky B, Neuert G, Oudenaarden A (2012) Using gene expression noise to understand gene regulation. Science 336: 183–187CrossrefCASPubMedWeb of Science®Google Scholar
- Narula J, Tiwari A, Igoshin OA (2016) Role of autoregulation and relative synthesis of operon partners in alternative sigma factor networks. PLoS Comput Biol 12: e1005267CrossrefPubMedWeb of Science®Google Scholar
- Norman TM, Lord ND, Paulsson J, Losick R (2013) Memory and modularity in cell-fate decision making. Nature 503: 481–486CrossrefCASPubMedWeb of Science®Google Scholar
- Novick A, Weiner M (1957) Enzyme induction as an all-or-none phenomenon. Proc Natl Acad Sci USA 43: 553–566CrossrefCASPubMedWeb of Science®Google Scholar
- Nyström T (2004) Growth versus maintenance: a trade-off dictated by RNA polymerase availability and sigma factor competition? Mol Microbiol 54: 855–862Wiley Online LibraryCASPubMedWeb of Science®Google Scholar
- Park J, Dies M, Lin Y, Hormoz S, Smith-Unna SE, Quinodoz S, Hernández-Jiménez MJ, Garcia-Ojalvo J, Locke JCW, Elowitz MB (2018) Molecular time sharing through dynamic pulsing in single cells. Cell Syst 6: 216–229.e15CrossrefCASPubMedWeb of Science®Google Scholar
- Patange O, Schwall C, Jones M, Villava C, Griffith DA, Phillips A, Locke JCW (2018) Escherichia coli can survive stress by noisy growth modulation. Nat Commun 9: 5333CrossrefCASPubMedWeb of Science®Google Scholar
- Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C (2017) Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods 14: 417–419CrossrefCASPubMedWeb of Science®Google Scholar
- Pinto D, Vecchione S, Wu H, Mauri M, Mascher T, Fritz G (2018) Engineering orthogonal synthetic timer circuits based on extracytoplasmic function σ factors. Nucleic Acids Res 46: 7450–7464CrossrefCASPubMedWeb of Science®Google Scholar
- Rackauckas C, Nie Q (2017a) Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling with memory. Discrete Continuous Dyn Syst Ser B 22: 2731–2761CrossrefPubMedWeb of Science®Google Scholar
- Rackauckas C, Nie Q (2017b) DifferentialEquations.jl – a performant and feature-rich ecosystem for solving differential equations in Julia. J Open Res Softw 5: 15CrossrefGoogle Scholar
- Raj A, Oudenaarden A (2008) Nature, nurture, or chance: stochastic gene expression and its consequences. Cell 135: 216–226CrossrefCASPubMedWeb of Science®Google Scholar
- Rhodius VA, Segall-Shapiro TH, Sharon BD, Ghodasara A, Orlova E, Tabakh H, Burkhardt DH, Clancy K, Peterson TC, Gross CA et al (2013) Design of orthogonal genetic switches based on a crosstalk map of σs, anti-σs, and promoters. Mol Syst Biol 9: 702Wiley Online LibraryCASPubMedWeb of Science®Google Scholar
- Rotem E, Loinger A, Ronin I, Levin-Reisman I, Gabay C, Shoresh N, Biham O, Balaban NQ (2010) Regulation of phenotypic variability by a threshold-based mechanism underlies bacterial persistence. Proc Natl Acad Sci USA 107: 12541–12546CrossrefCASPubMedWeb of Science®Google Scholar
- Runde S, Molière N, Heinz A, Maisonneuve E, Janczikowski A, Elsholz AKW, Gerth U, Hecker M, Turgay K (2014) The role of thiol oxidative stress response in heat-induced protein aggregate formation during thermotolerance in Bacillus subtilis. Mol Microbiol 91: 1036–1052Wiley Online LibraryCASPubMedWeb of Science®Google Scholar
- Soneson C, Love MI, Robinson MD (2015) Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4: 152CrossrefPubMedGoogle Scholar
- Spizizen J (1958) Transformation of biochemically deficient strains of Bacillus subtilis by deoxyribonucleate. Proc Natl Acad Sci USA 44: 1072–1078CrossrefCASPubMedWeb of Science®Google Scholar
- Sureka K, Ghosh B, Dasgupta A, Basu J, Kundu M, Bose I (2008) Positive feedback and noise activate the stringent response regulator rel in mycobacteria. PLoS One 3: e1771CrossrefCASPubMedWeb of Science®Google Scholar
- Taheri-Araghi S, Bradde S, Sauls JT, Hill NS, Levin PA, Paulsson J, Vergassola M, Jun S (2015) Cell-size control and homeostasis in bacteria. Curr Biol 25: 385–391CrossrefCASPubMedWeb of Science®Google Scholar
- Tiwari A, Ray JCJ, Narula J, Igoshin OA (2011) Bistable responses in bacterial genetic networks: designs and dynamical consequences. Math Biosci 231: 76–89CrossrefPubMedWeb of Science®Google Scholar
- Tjalsma H, Bolhuis A, Roosmalen ML, Wiegert T, Schumann W, Broekhuizen CP, Quax WJ, Venema G, Bron S, Dijl JM (1998) Functional analysis of the secretory precursor processing machinery of Bacillus subtilis: identification of a eubacterial homolog of archaeal and eukaryotic signal peptidases. Genes Dev 12: 2318–2331CrossrefCASPubMedWeb of Science®Google Scholar
- Veening J-W, Smits WK, Kuipers OP (2008a) Bistability, epigenetics, and bet-hedging in bacteria. Annu Rev Microbiol 62: 193–210CrossrefCASPubMedWeb of Science®Google Scholar
- Veening J-W, Stewart EJ, Berngruber TW, Taddei F, Kuipers OP, Hamoen LW (2008b) Bet-hedging and epigenetic inheritance in bacterial cell development. Proc Natl Acad Sci USA 105: 4393–4398CrossrefCASPubMedWeb of Science®Google Scholar
- Wang P, Robert L, Pelletier J, Dang WL, Taddei F, Wright A, Jun S (2010) Robust growth of Escherichia coli. Curr Biol 20: 1099–1103CrossrefCASPubMedWeb of Science®Google Scholar
- Xiong W, Ferrell JE (2003) A positive-feedback-based bistable “memory module” that governs a cell fate decision. Nature 426: 460–465CrossrefCASPubMedWeb of Science®Google Scholar
- Young JW, Locke JCW, Altinok A, Rosenfeld N, Bacarian T, Swain PS, Mjolsness E, Elowitz MB (2011) Measuring single-cell gene expression dynamics in bacteria using fluorescence time-lapse microscopy. Nat Protoc 7: 80–88CrossrefCASPubMedWeb of Science®Google Scholar
- Young JW, Locke JCW, Elowitz MB (2013) Rate of environmental change determines stress response specificity. Proc Natl Acad Sci USA 110: 4140–4145CrossrefCASPubMedWeb of Science®Google Scholar