Transient heterogeneity in extracellular protease production by Bacillus subtilis.

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.

Introduction

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 ). Activation of the competence or sporulation pathway only occurs in part of the population (Smits ). Since there are clearly two distinguishable cell types in both cases, this phenotypic variation was described as exhibiting bistability (Fujita ; Maamar and Dubnau, 2005; Smits ; Veening ).

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 ; Mader ). 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 ). For DegU to activate aprE gene expression, it needs to be phosphorylated by the cytoplasmic DegS sensor protein (Mukai ; Dahl ). 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 ; Msadek, 1999). Interestingly, all of these regulators are under direct or indirect negative control of the key sporulation regulator Spo0A∼P (Fujita ). 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).

Figure 1

aprE activation by a logic AND gate. (A) Kinetic scheme showing the reactions used to generate the computational model. Dotted boxes show the DegU system, which includes positive feedback and the bistable sporulation switch. Thin arrows represent basal transcription reactions, whereas thick arrows represent activated levels. For simplicity of the scheme, several reactions are omitted (e.g. SinR binding/dissociation, degradation of mRNA/proteins). (B) Cells require both DegU∼P AND Spo0A∼P to activate aprE transcription.

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.

Results

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 ; Fujita ; Veening ). 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 ; Eichenberger ). 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 ). 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.

Figure 2

Stationary phase heterogeneity in B. subtilis. Strains acoA-gfp, hag-gfp, srfA-gfp, sunA-gfp, bpr-gfp and aprE-gfp were grown in TY medium and cells were collected for phase-contrast and fluorescence microscopy. Time is given in hours relative to the transition point between exponential and stationary growth phases (T0). Left panels are from T=−2, middle panels T=5 and right panels T=20. Images are an overlay of both channels. Arrows indicate endospore-containing cells.

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 (P) fused to gfp and the upstream region of degU (P), 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 ; 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.

Figure 3

Heterogeneous expression of degU. Strains (A) A:degS-gfp and (B) A:degSU-gfp were grown in TY medium and cells were collected for phase-contrast and fluorescence microscopy. Time is given in hours relative to the transition point between exponential and stationary growth phases (T0). Images are an overlay of both channels. Arrows indicate heterogeneously expressing cells. (C) Quantification of fluorescence microscopy images. Micrographs were quantified and plotted as histograms. Every bin represents a certain percentage of the total population with an average fluorescence level in arbitrary units (AU) within that bin (N=cell count).

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 P-gfp is 126 arbitrary fluorescence units (AU) above background levels compared to 675 AU for P-gfp. Furthermore, the exposure times required to capture P-gfp and P-gfp is twice compared to that used to capture P-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 ). 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 ; Dartois ). 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.

Figure 4

Positive feedback within the DegU system. Strains (A) degS-gfp/hy, (B) degU-gfp/hy/ΔdegU and (C) degU-gfp/hy were grown in TY medium and induced with 1 mM IPTG at the mid-exponential growth phase. After 3 h of induction, cells were collected for phase-contrast and fluorescence microscopy. Images are an overlay of both channels. Micrographs were quantified and plotted as histograms. Every bin represents a certain percentage of the total population with an average fluorescence level above background in arbitrary units (AU) within that bin (N=cell count).

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).

Figure 5

Activation of aprE by DegUhy. Strains (A) aprE-gfp/hy and (B) aprE-gfp/hy/ΔSU were grown in TY medium and induced with IPTG at the mid-exponential growth phase. After 2 h of induction, cells were collected for flow cytometric analyses. At least 12 000 cells were measured for each strain. (C) Mean fluorescence of three independent experiments is plotted for the indicated IPTG concentrations. Error bars indicate the standard error of the mean. Strain aprE-gfp/hy is indicated as circles (•) and aprE-gfp/hy/ΔSU is shown as open triangles (▵). (D) The percentage of cells with a fluorescence level higher than channel 530 is depicted as a function of IPTG concentration. This cutoff was chosen arbitrarily to highlight the differences of highly expressing cells between the strains. The average value of three independent experiments is plotted. Error bars indicate the standard error of the mean.

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 ; Kunst ; 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 ). Deactivation of SinR is achieved by Spo0A∼P-dependent activation of the antagonist SinI (Bai ). It was shown that aprE expression is strongly reduced in a spo0A mutant background (Ferrari , 1988; Olmos ). 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 ; Gonzalez-Pastor ), 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.

Figure 6

Intertwinement between the Spo0A and DegU circuitries. (A) Strains 168, aprE-gfp, aprE-gfp/Δ0A, aprE-gfp/ΔdegU and aprE-gfp/ΔsigF were grown overnight in TY medium and cells were analyzed by flow cytometry. (B) Strains aprE-gfp/hy and aprE-gfp/hy/ΔabrB were grown in TY medium and induced with 1 mM IPTG at the mid-exponential growth phase. Cells were collected for flow cytometric analyses after 1 and 2 h of induction. At least 12 000 cells were measured for each strain and representative flow cytometry profiles are shown.

Table 1

Mean fluorescence as determined by flow cytometry

Strain	Mean fluorescence (AU)a
168 (no gfp)	6.2±0.1
aprE-gfp	16.7±3.3
aprE-gfp, Δspo0A	6.7±0.2
aprE-gfp, ΔdegU	7.2±0.2
aprE-gfp, ΔsigF	12.2±1.8

aFluorescence is in arbitrary units (AU). The data are based on eight biological replicates from different cultures measured on different days. For every measurement, at least 9000 cells were analyzed. The range (±) indicates the standard error of the mean.

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 ; Veening ). 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 ), 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.

Table 2

Parameters used for the simulations

Parameter	Value	Description
kmdeg	10−2 s−1	Rate of mRNA degradation
kT	4 × 10−2 s−1	Rate of translation per mRNA
kpdeg	4 × 10−4 s−1	Rate of protein degradation
I0	4 × 10−3 #/s	Transcription rate of degU from inactivated promoter
Im	4.8 × 10−2 #/s	Transcription rate of degU from activated promoter
kd	0.1 s−1	Rate constant of DegU dimer dissociation
kU	Equation (8)	Rate of degU mRNA transcription
kdeph	0.05 s−1	Rate constant of DegU dephosphorylation
kph	6 × 10−3 s−1	Rate constant of DegU phosphorylation
K1	7 # (degU) 10 # (aprE)	Equilibrium constants of DegU dimer binding to DNA promoters
kE	Equation (9)	Rate of aprE mRNA transcription
Kr1	7 #	Dissociation constant of AbrB binding to DNA promoters
Kr2	7 #	Dissociation constant of SinR binding to DNA promoters
IE0	2 × 10−2 #/s	Transcription rate of aprE from inactivated promoter
IEm	0.4 #/s	Transcription rate of aprE from activated promoter

aWe assume a cell volume of 10−15 l and used copy numbers (#) as a unit for effective concentration. 602 # corresponds to 1 μM/l concentration.

Figure 7

The mechanism responsible for heterogeneous aprE expression as suggested by modeling. (A) Deterministic model predicts bistability in the level of DegU protein for the range of phosphorylation rates shown. Within that range, depending on initial conditions, the system can be in either of two stable steady states—high-DegU (DegU-ON) and low-DegU (DegU-OFF) states. For values above the boundary, the system can be stable only in the ON state. (B, E) Stochastic simulations demonstrating highly variable switching kinetics into stable ON state from initially OFF state for DegU (B) and AprE (E). The parameter kph was set to the value 0.006 s−1 indicated by * in panel A. Initially, the concentrations are at equilibrium for kph=0.001 s−1 (OFF state). The dashed line represents the result of deterministic simulations. For panels C–F, an ensemble of 5000 cells was simulated. (C) Distribution of DegU levels at different times. (D) Distribution of the repressors used in our simulation. (F) Distribution of aprE-expressing cells at different times after window opening. Note that the first peak in the histogram is not drawn to scale because most of the cells with AprE OFF are in the first bin and that this peak is more spread in actual in vivo single-cell measurements due to autofluorescence.

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 ). 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 ; Veening , 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.

Figure 8

Single-cell aprE expression trajectories. Individual cells (P-gfp) were followed by time-lapse microscopy. (A) Fluorescence of all cells within the microcolony at the indicated time (from the beginning of the time lapse; see Supplementary Movie S1). Every bin represents the percentage of the total population with an average fluorescence level above background in arbitrary units (AU) within that bin (N=cell count). (B–D) Three sibling pairs were followed in time, after their birth. Circles depict the ‘left' and triangles the ‘right' sister cell. Fluorescence is indicated in arbitrary units above background.

Discussion

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 ). 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 ). 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 ; Maamar ). 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 ). 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 ). 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 ). 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 ) 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 , 2008) and are explained in more detail in Supplementary information.

Subpopulation fractionation

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 ). 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 ) 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.

Transcriptome analysis

DNA-microarrays containing specific oligonucleotides for all 4107 open reading frames (ORF) of B. subtilis were prepared as described (Lulko ). RNA was extracted and prepared for hybridization as described (Lulko ). 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 ) (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 ). Microarray data are accessible from the publicly available Gene Expression Omnibus repository (http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE9266.

Modeling procedures

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 p 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 ) and MATLAB© (Mathworks Inc.). An SBML file of the model can be downloaded from the journal website.

Supplementary Figures S1–S7

Supplementary Movie S1

Supplementary Tables S1–S3