The most sophisticated survival strategy Bacillus subtilis employs is the differentiation of a subpopulation of cells into highly resistant endospores. To examine the expression patterns of non‐sporulating cells within heterogeneous populations, we used buoyant density centrifugation to separate vegetative cells from endospore‐containing cells and compared the transcriptome profiles of both subpopulations. This demonstrated the differential expression of various regulons. Subsequent single‐cell analyses using promoter‐gfp fusions confirmed our microarray results. Surprisingly, only part of the vegetative subpopulation highly and transiently expresses genes encoding the extracellular proteases Bpr (bacillopeptidase) and AprE (subtilisin), both of which are under the control of the DegU transcriptional regulator. As these proteases and their degradation products freely diffuse within the liquid growth medium, all cells within the clonal population are expected to benefit from their activities, suggesting that B. subtilis employs cooperative or even altruistic behavior. To unravel the mechanisms by which protease production heterogeneity within the non‐sporulating subpopulation is established, we performed a series of genetic experiments combined with mathematical modeling. Simulations with our model yield valuable insights into how population heterogeneity may arise by the relatively long and variable response times within the DegU autoactivating pathway.
When environmental conditions become unfavorable, the Gram‐positive model bacterium Bacillus subtilis is able to employ a number of adaptive responses, such as competence development for DNA uptake, endospore formation and the production and secretion of proteolytic enzymes (Dubnau and Lovett, 2002; Piggot and Losick, 2002; Tjalsma et al, 2004). Activation of the competence or sporulation pathway occurs only in part of the population (Cahn and Fox, 1968; Hadden and Nester, 1968; Chung et al, 1994). As there are clearly two distinguishable cell types in both cases, this phenotypic variation was described as exhibiting bistability (Fujita et al, 2005; Maamar and Dubnau, 2005; Smits et al, 2005; Veening et al, 2005).
As most gene expression experiments in the stationary growth phase have been performed on the basis of population‐wide studies, little is known about the gene expression profiles of the specific subpopulations. Previous genome‐wide studies on sporulating cultures might have masked non‐sporulation‐related gene expression and putative additional levels of heterogeneity (Fawcett et al, 2000; Eichenberger et al, 2004). To reveal the expression patterns of non‐sporulating cells within isogenic sporulating cultures of B. subtilis, both a genome‐wide and a single‐cell approach were used. First, we developed a method to separate endospore‐containing cells from vegetative cells using buoyant density gradient centrifugation. The transcriptomes of the resulting subpopulations were compared using DNA‐microarray technology. This analysis revealed the occurrence of substantial heterogeneity in gene expression patterns within the isogenic B. subtilis culture. Cells either sporulate or activate a number of adaptive regulatory networks such as motility and competence development. Subsequent single‐cell analyses using promoter‐GFP fusions confirmed the microarray results and, surprisingly, revealed further heterogeneity within the non‐sporulating subpopulation (Figure 2). Only part of the cells within the vegetative subpopulation highly and transiently expresses genes encoding the extracellular proteases Bpr (bacillopeptidase) and AprE (subtilisin), both of which are under the control of the DegU transcriptional regulator (Figure 1A). These extracellular proteases are known to act as scavenging enzymes and degrade large, complex proteins to smaller peptides, which can subsequently be taken up again as a new source of nutrients. Since these proteases and the products of their degradation activity freely disperse within the (liquid) growth medium, all cells within the clonal population are expected to benefit from their activities, indicating that B. subtilis might employ cooperative or altruistic behavior.
To obtain more mechanistic insights into how such heterogeneity is generated, we performed a series of genetic experiments and developed a quantitative mathematical model that accurately describes our (genetic) observations. Our data show that degU transcription is heterogeneous and gradually increases with time in a unimodal distribution. Autoactivation of DegU is critical in reaching high levels of aprE transcription. Furthermore, phosphorylation of Spo0A, the master sporulation transcription factor, is also required to relieve the aprE promoter from repression of a number of transcriptional regulators (e.g. AbrB and SinR). Our experimental data suggest a time‐window model of how heterogeneity of aprE gene expression is generated. The temporal window for aprE expression opens when at least two conditions are satisfied. Firstly, environmental signals result in phosphorylation of Spo0A and de‐repression of negative repressors such as AbrB and SinR. Because only part of the population reaches the Spo0A∼P levels that are required to relieve the aprE promoter, only part of the population is primed to activate aprE gene expression. Secondly, an increase in the DegU phosphorylation rate (or decrease in the dephosphorylation rate) results in a higher probability of activating the degU autostimulatory loop. Thus, only cells that have high levels of both Spo0A∼P and DegU∼P will highly express aprE (Figure 1B). The aprE expression window closes when the gradual increase of Spo0A∼P reaches the level required to initiate endospore formation (Fujita and Losick, 2005) and may also be affected by additional factors such as cell death or induced DegU proteolysis.
Using this information, we built a qualitative mathematical model that constitutes a logic AND circuit involving the bistable sporulation pathway and the DegU autoactivation pathway (DegU system). Our model and experimental data support the hypothesis that, in the late stationary phase, the DegU system functions with parameters where only the fully activated state is operational. However, positive feedback in the system with the potential to demonstrate bistability results in slow and stochastically variable transitions from the ‘OFF’ to the ‘ON’ state. Simulations of the integrative model yield valuable insights into how aprE population heterogeneity arises from the relatively long and variable response times within the DegU system and makes testable experimental predictions.
Within a late stationary phase culture of Bacillus subtilis, three distinct subpopulations of cells can be distinguished: 1) cells that sporulate, 2) cells that do not, or only to low levels express the extracellular proteases Bpr and AprE, 3) cells that highly express extracellular proteases.
Heterogeneous expression of these proteases is established by a logic AND gate via the DegU and Spo0A regulatory proteins.
A mathematical model was constructed which accurately describes our genetic observations and makes testable predictions of the system.
Advanced time‐lapse microscopy confirms our modeling predictions on variable and long response times in aprE activation.
When the Gram‐positive Bacillus subtilis bacterium encounters nutrient limitations and enters the stationary growth phase, it can engage in multiple adaptation strategies such as the secretion of extracellular proteases, biofilm formation, uptake of exogenous DNA (competence) or production of highly resistant endospores (Dubnau and Lovett, 2002; Piggot and Losick, 2002; Tjalsma et al, 2004). Activation of the competence or sporulation pathway only occurs in part of the population (Smits et al, 2006). Since there are clearly two distinguishable cell types in both cases, this phenotypic variation was described as exhibiting bistability (Fujita et al, 2005; Maamar and Dubnau, 2005; Smits et al, 2005; Veening et al, 2005).
A key feature of B. subtilis is its ability to produce large quantities of extracellular proteases (exoproteases). Both aprE (subtilisin) and bpr (bacillopeptidase) genes encode such proteolytic enzymes that are secreted and degrade proteins from the environment (Msadek, 1999). Besides proteins present in the growth medium, these proteases can also act on fragments released from (dead) cells present within the stationary phase culture. Expression of both aprE and bpr is under the control of the DegS–DegU two‐component system (Ogura et al, 2001; Mader et al, 2002). Recent work has shown that DegU is a master regulator for multicellular behavior and that there are categories of genes within its regulon that respond to either high or low levels of DegU∼P and the genes coding for exoproteases require high levels of DegU∼P (Kobayashi, 2007; Verhamme et al, 2007). For DegU to activate aprE gene expression, it needs to be phosphorylated by the cytoplasmic DegS sensor protein (Mukai et al, 1990; Dahl et al, 1991). The exact biochemical nature of the signal that causes activation of DegS remains elusive, although it has been shown that the regulatory cascade can be activated under conditions of salt stress (Kunst and Rapoport, 1995). In vivo results indicate that DegU is capable of forming homo‐dimers (Sanchez and Olmos, 2004). Sequence and transcription analyses indicate that there are two putative binding sites for the DegU dimer within the aprE promoter (http://dbtbs.hgc.jp/; Shimane and Ogura, 2004). In addition to positive regulation by DegU∼P, aprE is under direct negative control of at least three other transcriptional regulators (AbrB, SinR and ScoC) (Msadek et al, 1993; Msadek, 1999). Interestingly, all of these regulators are under direct or indirect negative control of the key sporulation regulator Spo0A∼P (Fujita et al, 2005). Thus, when nutrients become limiting, signals that feed into the sporulation pathway will, besides regulating initiation of spore formation, also result in de‐repression of the aprE promoter. Taken together, it seems that the aprE gene‐regulatory network can be envisaged as a logic AND circuit, for which a threshold level of both dimerized DegU∼P and Spo0A∼P is required to activate gene expression (Figure 1B).
Presence of such an intriguing network motif led us to investigate the expression profile of the aprE‐activating system. Using both genome‐wide and single‐cell analyses, we show that within a liquid stationary phase culture of isogenic B. subtilis cells, at least three different cell types can be distinguished: (1) endospore‐containing cells, (2) cells with no or low levels of aprE/bpr expression (‘OFF’) and (3) cells with high levels of aprE/bpr expression (‘ON’). To understand the mechanisms and dynamics that underlie heterogeneity in aprE expression, we have built an integrative mathematical model based on our new results combined with previously published results. We demonstrate that the observed stationary phase multistability (‘ON’, ‘OFF’ and endospore formers) can be generated by the underlying heterogeneity in the levels of the two master regulators of protease production: Spo0A∼P and DegU∼P. This indicates that B. subtilis can utilize a remarkably simple strategy to generate multistability without the necessity to create a complex switch with multiple steady states. Furthermore, our mathematical study demonstrates how a transcriptional regulator can govern a time‐dependent heterogeneous expression of its target genes.
Identifying stationary phase heterogeneity
Previously, it has been shown that only part of a B. subtilis culture forms endospores and that the positive feedback architecture of the sporulation regulatory pathway is pivotal in this bistable development (Chung et al, 1994; Fujita et al, 2005; Veening et al, 2005). As most high‐throughput gene expression experiments in the stationary growth phase have been performed on the basis of population‐wide studies, little is known about the gene expression profiles of the specific subpopulations (Fawcett et al, 2000; Eichenberger et al, 2004). To investigate stationary phase heterogeneity of B. subtilis, we developed a method using buoyant density gradient centrifugation to separate endospore‐containing cells from vegetative cells (see Materials and methods). Total RNA was extracted from both fractions (vegetative cells and endospore‐containing cells) and used in subsequent DNA‐microarray experiments. Genes that were significantly overexpressed in one fraction compared to the other were analyzed using the FIVA software (Blom et al, 2007). This particular software aids in identifying functional profiles based on gene expression data and generates an overview of affected biological processes showing which gene‐regulatory pathways are differentially expressed (Supplementary Figure S1). The endospore fraction differentially expresses only a small number of genes compared to vegetative cells, which could still be detected by our microarray experiment (see Materials and methods). The vegetative cells however appear to express several different gene‐regulatory pathways. First of all, the Spo0A (initiation of sporulation) and SigH and SigE (early sporulation sigma factors) regulons are upregulated in the average of the vegetative cell fraction compared to the endospore fraction. This suggests that parts of the vegetative cells are still in the process of forming an endospore. Secondly, vegetative cells within the sporulating culture have activated members of the SigD (motility) and ComA (competence development and surfactin production) regulons. Non‐spore formers were also shown to express genes involved in the utilization of overflow metabolites (CcpA, SigL regulon). However, the most significant regulon that was overrepresented in the vegetative subpopulation was the DegU regulon (Supplementary Figure S1).
To validate the DNA‐microarray experiments, some of the top hits from each category were tested for their putative heterogeneous expression. Therefore, promoter‐GFP fusions, which allow single‐cell analyses, were constructed (Figure 2). As expected, we observed heterogeneous expression of the genes acoA (acetoin uptake, CcpA/SigL regulon), sunA (lantibiotic, AbrB/ComA/Rok regulon), srfA (surfactin synthetase and competence development, ComA regulon) and hag (flagellin protein, SigD regulon). In all of these cases, expression was either not present or reduced in endospore‐containing cells after overnight growth (Figure 2, arrows). Expression from the acoA promoter was activated when cells entered stationary phase and this occurred in all cells of the population. The sunA and srfA reporters, however, were already activated during the exponential growth phase. The hag reporter was expressed in a bistable fashion during the exponential growth phase, and chained exponentially growing cells do not express hag whereas single cells do, as reported before (Kearns and Losick, 2005). When cells entered the stationary growth phase, however, all cells express GFP from the hag promoter.
Reporter strains of two members of the DegU regulon, bpr and aprE, stood out compared to the other strains. Similar to the other strains, expression of GFP from the aprE and bpr promoters was not, or only to low levels, present in endospore‐containing cells. However, in contrast to the other reporter strains, only part of the vegetative subpopulation highly expresses bpr and aprE (Figures 2 and 3C). Thus, in stationary phase cultures, three distinct cell types could be identified: (1) highly expressing aprE or bpr cells, 2) cells that do not (or only at low levels) express aprE or bpr and (3) sporulating cells (that do not, or not anymore, express aprE or bpr) (Figure 2). The ratio of cells highly expressing aprE or bpr appeared to depend on the sporulation efficiency of the culture (see Results).
Heterogeneous expression of degU
To investigate whether the observed heterogeneity in aprE expression could originate from heterogeneous expression of its activator, promoter GFP reporter strains were constructed. Since degU lies downstream in an operon with degS, two strains were constructed: the upstream region of degS (PdegS) fused to gfp and the upstream region of degU (PdegU), including degS and PdegS, fused to gfp. As shown in Figure 3A, degS is homogeneously expressed and mainly activated in the stationary growth phase. degU is not, or only to low levels, expressed in exponentially growing cells. Interestingly, expression of degU is more heterogeneous compared to degS expression during the stationary growth phase (Figure 3B, arrows). Total DegU protein levels and degU RNA levels were shown to be highest at this point in growth (Hata et al, 2001; Kobayashi, 2007), consistent with our GFP expression data. While degS shows reduced expression after overnight growth in sporulating cultures, part of the population still strongly expresses degU.
Quantification of the micrographs clearly shows that degU is activated in the stationary growth phase and shows a unimodal but long‐tailed distribution (Figure 3C). Expression levels of degS and degU are significantly lower than aprE and are just above background levels. degU levels are significantly higher than degS (>2 times). As shown in Figure 3C, the highest bin of fluorescence of PdegU‐gfp is 126 arbitrary fluorescence units (AU) above background levels compared to 675 AU for PaprE‐gfp. Furthermore, the exposure times required to capture PdegU‐gfp and PdegS‐gfp is twice compared to that used to capture PaprE‐gfp cells. Importantly, the distributions of cells that express GFP from either degU or aprE are quite different. This suggests that besides activation of degU, another condition needs to be met to highly activate aprE expression.
Positive feedback within the DegU system
Using a genome‐wide approach, it was shown that overproduction of DegUhy, a constitutively active mutant protein of DegU, leads to induced degU (Mader et al, 2002). This suggests the presence of a positive feedback loop in degU activation. In fact, a putative DegU∼P binding site for degU is located downstream of the degS promoter (Msadek et al, 1990; Dartois et al, 1998). To test whether DegUhy overexpression was able to activate GFP production from the degS and/or degU promoter, we constructed a strain where degUhy is under the control of the isopropyl‐β‐d‐thiogalactosidase (IPTG)‐inducible hyper_spank promoter and introduced this construct in the relevant reporter strains. As shown by fluorescence microscopy, artificial induction of DegUhy was unable to activate transcription from the degS promoter and might even repress degS transcription (Figure 4A). In the absence of a wild‐type copy of degU, induction of DegUhy activates transcription from the degU promoter and all cells are homogeneously induced. When DegUhy is induced in the presence of degU, thus in the presence of positive feedback, a clear bistable expression pattern is observed in that part of the population does not express GFP and part of the population highly expresses GFP from the degU promoter. Taken together, these results indicate (1) the presence of positive feedback within the DegU system (Figure 4B) and (2) DegU∼P autostimulation leads to heterogeneous degU expression (Figure 4C). Thus, DegU∼P autostimulation might be important for heterogeneous aprE expression.
Maximal aprE expression requires DegU autostimulation
To examine how DegU autostimulation contributes to aprE heterogeneity, we introduced the IPTG‐inducible DegUhy construct in the aprE‐gfp reporter strain. The resulting strain (aprE‐gfp/hy) was grown in the presence of different concentrations of IPTG and examined by flow cytometry (Figure 5). Flow cytometry allows us to conveniently probe a large number of cells. The experiment was performed during the exponential growth phase to prevent putative influences derived from stationary phase population heterogeneity. As shown in Figure 5A, in the wild‐type background (with autostimulation), induction of DegUhy with IPTG activates aprE gene expression, and already at low inducer concentrations, part of the population highly expresses GFP. To examine whether degU autostimulation is crucial for aprE expression, we removed the positive feedback by replacing the native degSU locus with an antibiotic resistance marker. As shown in Figure 5B and C, aprE gene expression is also efficiently activated in the degSU mutant background, and the average fluorescence is only slightly less in the mutant background, especially at low levels of DegUhy induction. However, when we specifically look at the subpopulation of cells that highly express aprE‐gfp (>channel 530), the differences between wild type and the degSU mutant are significant at lower IPTG levels (P<0.05, t‐test, 25, 50, 75 and 100 μM IPTG) (Figure 5D). At longer induction times, however, aprE‐gfp expression levels are significantly higher in the degSU mutant (P<0.001, Student's t‐test; Supplementary Figure S5).
aprE heterogeneity is modulated by Spo0A∼P
If heterogeneity of degU expression is solely responsible for heterogeneity in aprE expression, why are the population distributions of the two reporter strain so different (Figure 3C)? Besides the positive regulation of aprE gene expression by DegU∼P, it was shown that aprE is negatively regulated by a number of transcriptional repressors (AbrB, ScoC and SinR). All of these are under direct or indirect negative control of Spo0A∼P, the master sporulation regulator (Msadek et al, 1993; Kunst et al, 1994; Msadek, 1999). Expression of abrB is repressed at low levels of Spo0A∼P and scoC in its turn is activated by AbrB (Perego and Hoch, 1988; Fujita et al, 2005). Deactivation of SinR is achieved by Spo0A∼P‐dependent activation of the antagonist SinI (Bai et al, 1993). It was shown that aprE expression is strongly reduced in a spo0A mutant background (Ferrari et al, 1986, 1988; Olmos et al, 1996). This indicates that in the natural situation, the presence of DegU∼P alone is not sufficient to activate aprE gene transcription.
Previous studies on aprE and degU regulation were all performed on whole populations and never at the single‐cell level. Since initiation of sporulation by Spo0A∼P only occurs in a subpopulation of cells (Chung et al, 1994; Gonzalez‐Pastor et al, 2003), we examined to which extent Spo0A contributes to the heterogeneous expression pattern of aprE. Therefore, we introduced a spo0A mutation into our single‐cell reporter strain (aprE‐gfp/Δ0A). Flow cytometry experiments demonstrate that without a functional spo0A gene, aprE is not or only to low levels expressed after overnight growth in TY medium (tryptone/yeast extract) (Figure 6A). As expected, a culture of a strain mutant for degU does not exhibit highly expressing aprE cells as well (Figure 6A). Actual spore formation is not required for aprE expression, as a strain mutant for SigF (aprE‐gfp/ΔsigF), a sporulation‐specific sigma factor that acts downstream of Spo0A, still significantly expresses aprE (Figure 6A). To quantify the differences between strains, we measured the average fluorescence of each strain of eight biological replicates using flow cytometry (Table I). Taken together, these results show that aprE‐gfp expression levels are significantly reduced in the spo0A and degU mutants.
To corroborate the intertwinement of the DegU with the Spo0A system, we tested the aprE‐gfp response to DegUhy in the absence of the negative regulator AbrB. As shown in Figure 6B, expression of GFP from the aprE promoter is more pronounced and more sensitive to induction of DegUhy in the abrB mutant background compared to the wild‐type strain. Induction of DegUhy in a spo0A mutant background was able to activate aprE gene expression although the response was delayed and reduced (Supplementary Figure S2). This shows that for high aprE expression, the Spo0A‐regulated repressors need to be de‐repressed.
Mathematical modeling predicts transient heterogeneous aprE expression
The simplest interpretation of our results is that heterogeneity in aprE expression is controlled by the intertwinement of two signal transduction cascades. First of all, to activate aprE, the negative control elements (i.e. AbrB, SinR and ScoC) need to be removed. This happens via activation of the sporulation regulator, Spo0A, which is under the control of a bistable switch (Fujita et al, 2005; Veening et al, 2005). Secondly, to activate aprE, cells need to activate degU via autostimulated transcription and for this DegU needs to be phosphorylated by DegS and dimerize (Figure 1A). We assume that a basal level of degU expression is present even in the non‐activated promoter (see Figure 3) and that the positive feedback of degU is controlled by DegS‐modulated phosphorylation of DegU (Supplementary Figure S5). Thus, cells that have highly activated aprE have reached a threshold level of both DegU∼P and Spo0A∼P. However, when Spo0A∼P levels continue to rise and get too high, cells commit to sporulation and aprE expression diminishes (see Discussion). It was shown that Spo0A∼P indirectly activates degU at low levels and represses it at high Spo0A∼P levels (Fujita et al, 2005), indicating the presence of more complex regulation than described here. For simplicity, these interactions were not taken into account. Simulations indicate that the putative coupling between degU expression and Spo0A∼P levels does not qualitatively alter the conclusions of the model (data not shown). Overall, the system seems to consists of two intertwining networks, which are both required to generate heterogeneity in aprE expression. Taking these parameters into account, we wondered whether an integrative mathematical model could predict the observed multistable phenotype (aprE ‘ON’, ‘OFF’ and endospores) on the basis of the proposed logic AND‐gate system. Specifically, we set out to uncover the mechanism by which the autoregulatory feedback loop results in a highly variable degU expression pattern, and how this interplays with the heterogeneity of Spo0A∼P‐controlled repressors AbrB and SinR to generate transiently heterogeneous production of exoproteases.
The essential reaction scheme used for modeling is shown in Figure 1A. Complete details of reactions, kinetics and the parameters used are summarized in Materials and methods and in Table II. The choice of parameter values and their ranges that result in bistable or monostable behavior of the DegU system are discussed in Supplementary information. It should be mentioned that very few of the parameters have been measured experimentally. Therefore, we concentrated on qualitative rather than quantitative predictions of the model. We first modeled the DegU autostimulatory network. We assume that the rate of increase in the DegU phosphorylation to dephosphorylation ratio (kph/kdeph) serves as a signal to activate degU transcription in stationary phase cells. This hypothesis is consistent with the data presented above. A deterministic (concentration based) simulation predicts a hysteretic bistable steady‐state response of DegU production that depends on the rate of DegU phosphorylation (Figure 7A). This means that for a range of dephosphorylation rates, two steady states are possible. The ‘DegU OFF’ state corresponds to inactivated levels of degU expression and the ‘DegU ON’ state corresponds to an almost maximal activation of the positive feedback loop. To predict the distribution of DegU in cells, we constructed a stochastic model of the DegU switch. Since very little is known about the biochemical parameter values characterizing the DegU system, we performed stochastic simulations for ranges of parameters that are in the same order as in other two‐component systems (see Supplementary information for more discussion). Because our experimental results (Figure 3) indicate a unimodal rather than a bimodal distribution of DegU protein in the wild‐type strains under natural conditions, we find that the model most reliably explains the experimental observations for the parameter values outside the bistability range but not too far from its boundary (e.g., as indicated by * in Figure 7A and Supplementary Figure S7). For these parameter values, the transitions from the OFF to ON state are more likely to occur than the opposite transition. Therefore, given enough time, all the cells with a phosphorylation rate above a certain threshold will activate the autostimulation. However, the time needed for such activation can be relatively long and highly variable as indicated by several sample trajectories from the stochastic simulations (Figure 7B). Therefore, we predict that stochastic variability in the time to activate the positive feedback loop in degU expression is the main source of DegU heterogeneity after overnight grown cultures. Figure 7C shows that the model predicts a highly variable DegU distribution within the population. We do not see a clear bimodal distribution in agreement with the experimentally measured distribution (Figure 3). The fraction of cells that highly express DegU increases with time.
Since a mechanistic model that describes Spo0A∼P bistability does not exist, we decided to use experimentally measured distributions. Previously, we have described a double‐labeled reporter strain that is indicative of AbrB and Spo0A∼P concentrations in single cells (Veening et al, 2004). In this strain, the gene encoding the yellow fluorescent protein (YFP) is under the control of the abrB promoter and the cyan fluorescent protein (CFP) is under the control of the spoIIA promoter. Activation of spoIIA is a good indication for cells that have accumulated Spo0A∼P (Fujita et al, 2005; Veening et al, 2005, 2008). Quantification of cells grown overnight in TY medium shows that approximately 43% of cells are activated in spoIIA transcription and not express abrB anymore, whereas approximately 34% of the population still expresses abrB and not yet expresses spoIIA. In these cultures, approximately 22% of cells have formed endospores and do not visibly express either CFP or YFP. Since the model does not include endospore formation, we excluded these events from the statistics and used a probability of approximately 45% (34/(34+43)∼45%) that the repressors are active in a cell. Using this information, we modeled a bimodal distribution of AbrB and SinR concentrations in single cells (Figure 7D). For simplicity, we limit the model to account for only two repressors (SinR and AbrB). Because activity of both repressors is under the control of Spo0A∼P, we assumed that the active repressor concentration varies proportionally to one another.
An integrated model where DegU heterogeneity is combined with the bimodal distribution of the repressors is successful in generating a heterogeneous aprE expression pattern. Sample trajectories from the stochastic simulations show that aprE expression is highly variable between individual cells (Figure 7E). Importantly, this model predicts that with increasing time, a larger population of the cells highly expresses aprE (Figure 7F).
Experimental validation of the model
To test our predictions in vivo, we monitored the single‐cell aprE gene expression pattern over an extended period of time. As predicted by the model, over time, a larger fraction of cells highly express aprE. However, after approximately 17 h of growth in vivo in liquid cultures, the aprE subpopulation does not increase further and rather decreases (Supplementary Figure S4). This is likely caused by cell death and spore formation within the culture. Since these events are not included in the model, the population of aprE‐expressing cells continues to increase with time in our simulations.
Since our model predicts long response times of degU activation, and subsequently aprE activation, we setout to determine single‐cell trajectories of degU expression. However, because of the low expression levels of degU, long exposure times are required, which induces photo‐toxicity and bleaching. This makes single‐cell time‐lapse microscopy with the degU reporter at this moment technically challenging. Instead, we focused on aprE expression and generated time‐lapse movies of the aprE‐gfp reporter strain where we follow single cells through time and division. Movies were made where single cells were grown into sporulating microcolonies and fluorescence was followed at timely intervals (see Supplementary Movie S1). Cells do not express aprE during the exponential growth phase, consistent with the results obtained in liquid media. When cells enter the stationary growth phase, only a subpopulation of cells activates aprE. With increasing time, more cells activate aprE (Figure 8A). In some cases, cells switch off expression of aprE after having the gene switched on (Supplementary Movie S1). Three sister pairs were followed from their birth and their individual GFP concentrations were plotted against time. Comparing siblings excludes epigenetic, temporal and geographic factors involved in the decision making to activate aprE. As shown in Figure 8, cells activate aprE in a highly variable fashion, consistent with our mathematical model.
Heterogeneity in stationary phase gene expression
Examining stationary phase B. subtilis cultures by both genome‐wide and single‐cell analyses revealed substantial heterogeneity. Genome‐wide transcriptome analyses of subpopulations separated on the basis of their buoyant density showed that non‐spore formers are either in the stage of becoming spore formers and/or have activated a large variety of members of other stationary phase processes including genes required for competence development, motility and the production of exoproteases. For at least srfA and aprE, it has been shown biochemically that they are transcribed in an σA‐dependent fashion and many more σA‐dependent genes were picked up in the microarray experiment. It has long been known that although the levels of this major housekeeping RNA polymerase sigma factor, σA, remains constant throughout growth and sporulation, the transcription levels of genes under its control decrease significantly upon entry into sporulation (Linn et al, 1973). It is unclear whether this decrease is due to σA activity regulation or due to its competition with alternative, sporulation‐specific sigma factors (Hicks and Grossman, 1996; Fujita, 2000). In any case, reduced σA activity of genes that are normally expressed during growth may very well be responsible for some of the heterogeneous expression patterns we have observed in stationary phase cultures, as only part of a culture sporulates.
One of the most interesting heterogeneously expressed genes is aprE, which encodes subtilisin, an (industrially) important exoprotease. After overnight growth, three distinct cell types could be identified: highly expressing aprE cells, cells that do not (or to a low level) express aprE and sporulating cells (that do not, or not anymore, express aprE). Flow cytometry analysis indicates a long‐tailed distribution of aprE‐expressing cells (Figure 6). aprE is under direct control of the DegS/DegU two‐component system. Therefore, we wondered whether degU is also heterogeneously expressed. Indeed, degU expression is also specifically activated during the stationary growth phase and shows a unimodal but long‐tailed distribution. A heterogeneous expression pattern is most apparent in sporulating cultures after overnight growth, in that some cells highly express degU whereas others do to a much lesser extent (Figure 3). Using synthetic gene‐regulatory networks, it was found that stochasticity in gene expression (noise), amplified by positive feedback, can play a major role in phenotypic variability and can generate bistability (Kærn et al, 2005). Experiments using an artificially inducible variant of a mutant DegU protein that is constitutively active (DegUhy) show that degU is subject to bistable regulation (Figure 4). In the absence of a functional degU gene, all cells homogenously express GFP from the degU promoter upon DegUhy induction, indicating the presence of positive feedback. Intriguingly, in the presence of a wild‐type copy of degU (i.e. when the positive feedback loop is present), a clear bistable expression pattern is observed upon DegUhy induction, with only part of the cells (extremely) highly expressing GFP from the degU promoter whereas other cells not expressing (Figure 4C). When we examine the total DegU protein levels under these conditions using western blotting, we observe a similar situation in that a visible increase in total DegU protein levels compared to a degU mutant could not be detected (Supplementary Figure S3). However, there was substantial proteolytic cleavage. Prolonged induction of DegUhy also results in reduced aprE‐gfp expression in the wild‐type strain but not in the degSU mutant (Supplementary Figure S5). These results suggest that the production of unphosphorylated DegU directly or indirectly represses degU transcription or that there is an as yet unknown mechanism responsible for our observations. Interestingly, although our experiments clearly show the presence of DegU autoactivation, this does not result in a bistable expression pattern of degU in the wild‐type strain (Figure 3). We interpret this result by hypothesizing that the rates of phosphorylation–dephosphorylation have the values at which the system is monostable. Genetic perturbations, however, can make the system bistable (Figure 4). However, DegU autoactivation and potential bistability are important in establishing a heterogeneous aprE expression pattern (Figure 5).
The mechanism behind aprE heterogeneity
For competence development in B. subtilis, it has been shown that shutdown of basal expression of the master competence regulator ComK sets a time‐window for cells switching into the competent state (K‐state) (Leisner et al, 2007; Maamar et al, 2007). Our experimental data suggest a time‐window model of how heterogeneity of aprE gene expression is generated. The temporal window for aprE expression opens when two conditions are satisfied. Firstly, environmental signals result in phosphorylation of Spo0A and de‐repression of the negative repressors such as AbrB and SinR. Because only part of the population reaches the Spo0A∼P levels that are required to relieve the aprE promoter, only part of the population is primed to activate aprE gene expression. Secondly, an increase in the DegU phosphorylation rate (or decrease in the dephosphorylation rate) results in a higher probability of activating the degU autostimulatory loop. Thus, only cells that have high levels of both Spo0A∼P and DegU∼P will highly express aprE. The aprE expression window closes when the gradual increase of Spo0A∼P reaches the level required to initiate endospore formation (Fujita and Losick, 2005) and may also be affected by additional factors such as cell death, downregulation of degU expression or induced DegU proteolysis. Notably, existence of a bistable response and a gradual temporal increase in a regulator are not mutually exclusive but can represent two manifestations of the same phenomenon. Indeed, when the system is on the verge of the switching threshold, many ‘failed’ attempts are required before the positive feedback is eventually fully activated. Therefore, bistable systems with parameter values near the switching threshold show a slow response (see below).
A mathematical model that integrates the conditions required to activate aprE was indeed successful in generating a heterogeneous output (Figure 7). Simulations of this model show that the DegU system acts with a long response time and displays highly variable switching times. The model predicts that with increasing time, more cells should activate aprE, which is indeed the case (Figure 8 and Supplementary Figure S4).
Time‐series microscopic analysis of the degU reporter strain grown in batch cultures shows that degU indeed is heterogeneously expressed, as predicted by the model. Quantification of the microscopic images from strains aprE‐gfp and degSU‐gfp shows an increase of these reporters with increasing time, consistent with the flow cytometry data (Figure 3). To directly test our model, which predicts stochastically variable activation of aprE expression, automated time‐lapse analyses using the aprE‐gfp reporter were performed. This demonstrated that activation of aprE indeed is highly variable among cells (Figure 8 and Supplementary Movie S1). It should be noted that the growth conditions within microcolonies are different compared to growth in liquid media, but it nevertheless demonstrates transient aprE expression as observed in liquid cultures.
A question that arises from our simulations is why the kinetics of DegU activation is so slow and noisy. Our simulations indicate that this result is not very sensitive to parameter values as long as the dephosphorylation rate is chosen to be on the verge of the bistability region. The results shown in Figure 7B, C, E and F are performed with a dephosphorylation rate slightly above the ‘ON’ threshold. Slow kinetics of bistable switches have been indicated by earlier studies (Savageau, 2002). The intuitive reason for such a response within the DegU system is as follows. Before full‐scale feedback activation, intrinsic and extrinsic noise in degU transcription, translation, phosphorylation and dimerization reactions will lead to transient formation of the DegU∼P dimer. This formation may lead to positive feedback firing, but only very briefly. Indeed with a basal level of DegU synthesis (before full feedback activation), the concentration of phosphorylated dimer is lower than the DNA dissociation constant. After each firing, more DegU is produced, but only part of these molecules is phosphorylated and an even smaller part forms dimers. These increases in DegU∼P dimers only slightly elevate the probability of the next firing of the positive feedback. Thus, many consecutive non‐probable reaction events need to occur before the DegU∼P dimer concentration reaches the level where it gives rise to a high probability of binding its autoactivating site. This explains the relatively slow response time of the system. Stochastic fluctuations, for instance, in the DNA‐binding rate will speed up the transition, whereas other fluctuations may slow down the transition. As the kinetics are determined by a number of random events, the system is highly sensitive to fluctuations and this explains the high variability in switching of the DegU system. A similar mechanism may be responsible for bistability and gradual activation of Spo0A∼P and perhaps other gene‐regulatory networks with a heterogeneous outcome.
Several examples of heterogeneity in isogenic bacterial populations with bimodal gene expression distributions have been reported (see Dubnau and Losick, 2006; Smits et al, 2006). In B. subtilis, besides sporulation bistability, it was shown that the alternative sigma factor SigD, which regulates motility, is not expressed in early exponentially growing cells and that the switch from the non‐motile to the motile state is subject to bistable regulation (Kearns and Losick, 2005). Competence development in B. subtilis was also shown to be a (transient) bistable process (Suel et al, 2006). The simulations and experiments reported here indicate one of the possible mechanisms to generate such bimodality: very slow and highly variable kinetics can be observed when the system parameters change just outside the bistable region and only the ON state is stable. This mechanism will result in an increase of the ON cells fraction as a function of time. In the DegU system, negative feedback or additional inputs such as downregulation of degU and high level of Spo0A∼P or induced DegU proteolysis may terminate the switching thereby locking the distribution.
Evolutionary role of heterogeneous exoprotease production
AprE and Bpr are scavenging proteins that are secreted into the growth medium and degrade (large) proteins into smaller peptides, which can be taken up and used as an alternative nutrient source (Msadek, 1999). As we report, only part of the vegetative population highly expresses aprE/bpr. Because small peptides can presumably be taken up by the entire population, cells that produce and secrete AprE/Bpr not only help themselves, but also all the cells within the liquid growth medium. This could reflect a simple form of altruism: a behavior that decreases the fitness of the altruistic individual while benefiting others (West et al, 2006). Heterogeneity in aprE gene expression ensures that not all cells commence into the costly production of AprE, but all cells within the clonal population benefit from the activity of the exoprotease. Future research will tell whether this is true altruism where the actors suffer from their behavior and the whole population benefits, or whether it is cooperative and that both the actor and the recipient mutually benefit.
Materials and methods
Strains and growth conditions
Strains and plasmids are listed in Supplementary Tables S1 and S2 respectively. Oligonucleotides are listed in Supplementary Table S3. Plasmid and strain construction are described in Supplementary information. TY medium contained tryptone (1%), yeast extract (0.5%), NaCl (0.5%), MnCl2 (0.1 mM) and NaOH (2.8 mM). Sporulation medium was prepared as described before (Schaeffer et al, 1965) and contained dehydrated nutrient broth (0.8%), NaOH (0.5 mM), MgSO4 (1 mM), KCl (1 g/l), Ca(NO3)2 (1 mM), MnCl2 (0.01 mM) and FeSO4 (0.001 mM). Strains were grown at 37°C and continuously shaken at 245 r.p.m. Overnight growth is typically an incubation time of 18 h after inoculation directly from a −80°C glycerol stock. To generate Figures 2 and 3, cultures were ‘balanced’ by growing them to mid‐exponential phase and subsequently re‐diluting into fresh medium. This protocol was repeated two times before the actual time‐series microscopy experiment was started. Because sporulation commences later in balanced cultures than in typical ‘overnight’ cultures, later time points (e.g. T=20) are shown for these experiments. When required, medium for Escherichia coli was supplemented with ampicillin (100 μg/ml); media for B. subtilis were supplemented with chloramphenicol (5 μg/ml), tetracycline (6 μg/ml), spectinomycin (100 μg/ml), erythromycin (5 μg/ml) or kanamycin (10 μg/ml). When indicated, IPTG was added to the medium. Microscopic and flow cytometric techniques were essentially performed as described before (Veening et al, 2005, 2008) and are explained in more detail in Supplementary information.
To separate spore formers from non‐spore formers, we made use of differences in buoyant densities of these cell types. Buoyant density gradient centrifugation on B. subtilis cells was successfully adopted previously to separate competent from non‐competent cells and ‘heavy’ spores from ‘light’ spores (Cahn and Fox, 1968; Dean and Douthit, 1974). Both methods are based on the differences in buoyant density of cells in different developmental stages. Since early spore formers contain either a phase‐dark or phase‐bright endospore and vegetative cells do not, it was to be expected that the presence of an endospore alters cellular buoyant density. It was shown previously that vegetative cells could be separated by mature free spores using this technique (Dean and Douthit, 1974; Siccardi et al, 1975). Buoyant density gradient centrifugation can be carried out by centrifugation on a medium with a high density and low viscosity. Renografin (or urografin) is a compound that was shown to exhibit no toxic effects on vegetative cells of B. megaterium (Tamir and Gilvarg, 1966). As renografin is no longer available on the market, we used Ultravist, a solution with similar characteristics as renografin. To separate spore formers (cells containing endospores) from non‐spore formers (vegetative cells), B. subtilis cells (strain iyfp‐abrB‐icfp‐IIA‐amyE; Veening et al, 2004) were grown in sporulation medium. Two and a half hours after cells had entered stationary growth and the first appearance of cells with an endospore could be identified, 20 ml of culture was centrifuged for 1 min at 12 000 r.p.m. and the cell pellet was resuspended in 3 ml of a 10% Ultravist‐370 solution (Schering Nederland B.V., Weesp, The Netherlands). This suspension was loaded on top of an ultracentrifuge tube (Ultra Clear, 14 × 89 mm) containing a discontinuous layered gradient of 37% (bottom layer) and 33, 31.5, 30 and 28.5% (top layer) of Ultravist‐370 solution. The loaded gradient was centrifuged for 1 h at 25 000 r.p.m. under vacuum at 4°C in a Beckman Ultracentrifuge (Beckman Coulter Nederland B.V., Mijdrecht, The Netherlands) using rotor SW41. After centrifugation, cells were extracted from the top layers (vegetative cells) and the pellet fraction (endospores), washed and frozen in liquid nitrogen and stored at −80°C, and fractions were later used for RNA extraction. Successful separation of the different cell fractions was verified by microscopy.
DNA‐microarrays containing specific oligonucleotides for all 4107 open reading frames (ORF) of B. subtilis were prepared as described (Lulko et al, 2007). RNA was extracted and prepared for hybridization as described (Lulko et al, 2007). Two biological replicates were performed and the DNA‐microarray slides contained two probes of each ORF. Thus, maximally four data points for each gene were obtained. It should be noted that due to the timely centrifugation treatment, RNA quality isolated from the recovered cells was compromised, indicating that the microarray data obtained will not completely reflect the true mRNA state of both subpopulations. Also, only a relatively small amount of signal (transcripts) could be detected for the endospore‐containing fraction. Because of this, standard microarray data analyses could not be effectively used on these data sets. Therefore, we undertook the following approach: raw data were subjected to MicroPrep analyses (van Hijum et al, 2003) (without a Lowess normalization) and genes with a coefficient of variance lower than 60% were taken into account for further analyses. Genes that showed a more than two‐fold overrepresentation in the vegetative fraction compared to the endospore fraction were considered to be downregulated (Supplementary Figure S1, ‘down’). Genes with a ratio (endospore/vegetative cells) higher than 1.6 were considered to be upregulated in the endospore fraction (Supplementary Figure S1, ‘up’) and genes with ratios between 0.5 and 1.6 were not considered to be significantly different between both fractions (Supplementary Figure S1, ‘not’). The resulting list was subsequently analyzed using the FIVA software (Blom et al, 2007). Microarray data are accessible from the publicly available Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE9266.
The main ingredients of the model are briefly summarized below:
1. To decrease the number of unknown parameters, we chose to model protein synthesis as two elementary reaction steps, transcription and translation. The rate of translation for both proteins is modeled as a product of rate constant kT to respective mRNA concentrations:
2. Rate of mRNA degradation is assumed to be the same for both messages:
3. Rates of degradation and dilution due to cell growth were assumed to be the same for all the proteins and protein complexes:
4. DegU can be phosphorylated or dephosphorylated:
5. Phosphorylated DegU reversibly forms dimers:
6. We assume that promoter binding kinetics is fast and, therefore, in quasi‐steady state.
7. For simplicity, we limit the model to account for only two repressors: SinR and AbrB. Because activity of both the repressors is under the control of Spo0A∼P, we simplified the model by assuming that active repressor concentrations vary proportionally to one another.
8. Expression of degU is activated by concentration of Dim ([Dim]) as follows:
Here I0 is the flux from inactivated promoter, Im is the flux from activated promoter and K1 is an equilibrium constant of dimer binding to the promoter site.
9. Expression of aprE is modulated by the level of the repressors and the level of the activator, DegU∼P dimer. By analyzing the sequence of the aprE promoter, we conclude that it is capable of binding two DegU dimers and we assume that both of these are required for activation. Sequence analysis predicts that binding of AbrB and DegU dimers is mutually exclusive. On the other hand, the SinR binding site is more than 200 bp upstream of the transcriptional start site and, therefore, assumed to be independent of AbrB and DegU binding:
Here pR0 is the probability that SinR is not bound:
p10 is the probability that AbrB is not bound and no more than one DegU dimer is bound:
p1D is the probability that two DegU dimers are bound:
10. As the biochemical details causing heterogeneity in the levels of Spo0A∼P in cells are not known, we introduced this heterogeneity into the model by drawing concentration of repressors from a bimodal distribution estimated by analyzing the expression of the abrB and spoIIA promoters (see Results). Proportional concentrations for both repressors were assumed.
Parameters used for simulations are summarized in Table II. Our choice of parameter values is discussed in Supplementary information. Deterministic and stochastic simulations were carried out using Dizzy (Ramsey et al, 2005) and MATLAB© (Mathworks Inc.). An SBML file of the model can be downloaded from the journal website.
J‐WV was supported by Grant ABC‐5587 from the Netherlands Organization of Scientific Research, Technology Foundation (NWO‐STW), by a Ramsay Fellowship from the Royal Netherlands Academy of Arts and Sciences (KNAW) and by a grant from the Biotechnology and Biological Sciences Research Council awarded to J Errington. LWH was supported by a Wellcome Trust Research Career Development Fellowship. OAI was supported by a start‐up fund from Rice University. We thank David Rudner for the generous gift of plasmid pDR111, Tarek Msadek for strains and Teruo Tanaka for strains and DegU antibodies. We thank Esther de Jong for excellent technical assistance, various members of our labs for critically reading the manuscript and the anonymous referees for useful comments.
Supplementary Figures S1–S7
Supplementary Movie S1
Supplementary Tables S1–S3
This is an open access article under the terms of the Creative Commons Attribution‐NonCommercial‐NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non‐commercial and no modifications or adaptations are made.
- Copyright © 2008 EMBO and Nature Publishing Group