An experimentally supported model of the Bacillus subtilis global transcriptional regulatory network

Abstract

Organisms from all domains of life use gene regulation networks to control cell growth, identity, function, and responses to environmental challenges. Although accurate global regulatory models would provide critical evolutionary and functional insights, they remain incomplete, even for the best studied organisms. Efforts to build comprehensive networks are confounded by challenges including network scale, degree of connectivity, complexity of organism-environment interactions, and difficulty of estimating the activity of regulatory factors. Taking advantage of the large number of known regulatory interactions in Bacillus subtilis and two transcriptomics datasets (including one with 38 separate experiments collected specifically for this study), we use a new combination of network component analysis and model selection to simultaneously estimate transcription factor activities and learn a substantially expanded transcriptional regulatory network for this bacterium. In total, we predict 2,258 novel regulatory interactions and recall 74% of the previously known interactions. We obtained experimental support for 391 (out of 635 evaluated) novel regulatory edges (62% accuracy), thus significantly increasing our understanding of various cell processes, such as spore formation.

Keywords: Bacillus subtilis; network inference; sporulation; transcriptional networks.

Figure 1. General workflow for inferring the B. subtilis transcription network

Two transcriptomic data compendia were used, one collected specifically for this study (strain PY79) and one previously published for strain BSB1 (Nicolas et al, 2012). Transcription factor activities (TFA) were estimated independently for each dataset using interactions in the gold standard (GS) extracted primarily from SubtiWiki (step 1). Datasets, estimated TFA, and priors on network structure (from the GS) were used as inputs for prediction of regulatory interactions (step 2). Next, output networks (one for each strain/dataset) were merged into a combined network (inferred TRN) (step 3) and prediction accuracy was evaluated in follow‐up experiments.

Figure 2. Incorporating Transcription Factor Activities (TFA) in the network inference procedure

1. Partial Pearson correlation between mRNA transcription levels (PY79 dataset) was computed for each TF–target gene pair in the GS (top histogram). Partial Pearson correlation was also computed between the estimated activity of a TF and the transcription of its targets (bottom histogram).

2. The advantage of estimating TFA is illustrated for three regulators. Each point corresponds to the results of one microarray experiment, and TFA are estimated for each experimental condition. Top panel: A nonlinear correlation is observed between comK transcription and transcription of ComK targets, whereas a strong linear correlation is obtained between ComK activity and transcription of ComK targets. Middle panel: No correlation is observed between codY transcription and transcription of CodY targets. CodY activity is modulated by GTP and branched chain amino acids (BCAA). A negative correlation is observed between estimated CodY activity and transcription of CodY targets. Bottom panel: Spo0A activity is modulated by phosphorylation. A better correlation is observed between Spo0A activity and transcription of Spo0A targets than between spo0A transcription and transcription of Spo0A targets.

Figure 3. Performance of network inference methods when incorporating TFA

1. Precision–Recall plot of the confidence‐ranked interaction networks using CLR, Genie3, and the Inferelator (no priors). Solid lines show performance using TFA. Dashed lines show performance when no TFA are used (when raw expression values for TFs are used as predictors). The numbers superimposed on each curve indicate the area under the curve.

2. Performance of BBSR‐TFA (AUPR: area under precision recall curve) on the combined, BSB1 and PY79, networks in the presence of false prior information. 50% of the edges in the GS are used as true priors, and various amounts of random edges are added. Performance is evaluated on the leave‐out set of interactions. Each point represents the median of five random samples of 50% of the GS set.

3. Support from KO data for the models predicted by BBSR‐TFA, Genie3 (G3), CLR, and a consensus method (META) that rank combines the prediction of the three methods. Methods were used without and with TFA (TFA tag). The number on top of each bar indicates the proportion of evaluated interactions with KO support for the corresponding method. Left and right panels show the support for each method when all interactions (recovered priors and novel interactions) and only novel interactions are considered, respectively.

Figure 4. Experimental support from KO data to the models predicted by BBSR, Genie3, and CLR

For each regulon with KO data, we assessed the proportion of predicted targets supported by the KO data. Results are presented for BBSR, Genie3 (G3), CLR, and a consensus method (META) that rank combines the prediction of the three methods. Methods were used without and with TFA (TFA tag). The number in parentheses next to the regulon's name indicates the total number of differentially transcribed genes in the corresponding KO data. The number on top of each bar indicates the proportion of evaluated interactions (recovered priors and novel interactions) supported by the KO data. This number is omitted if there was no significant (P‐value ≥ 0.01) enrichment for differentially transcribed genes in the predicted targets

Figure 5. Modular organization of the inferred TF network

Cytoscape view of the modular architecture of the inferred network, restricted to σ factors (octagons) and other transcription factors (circles). The size of each node reflects the total number of predicted targets. Modules are labeled based on the functional annotation of their members. Green edges are known interactions; blue edges are novel interactions. Pie chart within each node indicates the proportion of known members (i.e. present in the GS network, green) and novel members (blue) for each regulon. The corresponding Cytoscape file is provided in Dataset EV5.

Figure 6. Functional analysis of spore polysaccharide synthesis genes

1. Genomic organization of the ytdA, yfnH, and spsI gene regions (ytdA is a paralog of yfnH and spsI). Putative gene functions are color‐coded. β score for each prediction is indicated in parenthesis.

2. Left: Fluorescence microscopy images for the YtdA‐GFP fusion in sporulating cells in the indicated mutant backgrounds. Except where indicated otherwise, images were collected for sporulating cells at hour 6 after suspension in Sterlini–Mandelstam medium at 37°C. Right: Possible binding site for σ^K in the ytdA promoter. The consensus binding sequence for σ^K is also indicated (M is A or C).

3. Spore coat localization of YtdA‐GFP and SpsI‐GFP is dependent on cotE and independent of cotXYZ.

4. Subcellular localization of YtcB‐YFP, SpsJ‐YFP, and SpsK‐CFP during sporulation.

5. Time course and dual labeling analysis of YfnH‐YFP and SpsM‐CFP during sporulation.

6. Spore adhesion to glass: A spsI deletion mutant and a spsI ytdA double mutant display strong adhesion.

7. Spore adhesion to hydrocarbons (hexadecane): A spsI deletion mutant and a spsI ytdA double mutant display strong adhesion. A spsI yfnH double mutant and a spsI yfnH ytdA triple mutant display intermediate adhesion. Error bars represent the standard deviation for three independent experiments.