Divergence of the SigB regulon and pathogenesis of the Bacillus cereus sensu lato group

SigB binding site model building and regulon predictions

We began with nine 150bp DNA sequences encompassing the SigB-dependent promoters identified by Van Schaik et al.[44], expanded to 166 sequences by phylogenetic footprinting[46]. Redundant sequences were removed from this collection to yield a final training set of 130 sequences. The training set was then used in parallel to develop separate HMM and PWM models for the SigB-dependent promoters in the Bc sensu lato group. Our models are virtually identical to previously derived models for this DNA binding site in Bc ATCC14579[44] and B. subtilis[47] (Figure1). Each model was used to scan all genomes for potential SigB binding sites, and this information was coupled with transcriptome-derived transcriptional unit (TU) predictions to arrive at a predicted regulon structure for each organism.

Horizontal gene transfer leading to insertion of novel genes, gene duplication and divergence, or gene deletion events, are thought to be the most common mechanisms responsible for remodeling bacterial genome function. Collectively, these events can be easily identified as the insertion or deletion (indel) of protein coding sequence, although the underlying mechanism responsible for the appearance of an indel is not so easily surmised. Surprisingly, preliminary manual comparison of the predicted SigB regulons of the Bc sensu lato group suggested that the appearance of indels was not the predominant mechanism of divergence of these SigB regulons. In the majority of instances, differences between SigB regulon structure appeared to arise from the placement of SigB-specific promoters in the 5’ region of TUs that were shared among these organisms. That is, SigB regulon divergence appeared to rely on a process of promoter ‘re-assignment’ of genes from a common gene pool into/out of the SigB regulon. This promoter re-assignment appears to predominate over instances of indels that changed the total coding capacity of a given SigB regulon. To get a better quantitative measure of this impression, we performed pairwise comparisons of the predicted SigB regulons of all Bc sensu lato genomes (Figure2), by calculating a simple ratio of the number of predicted promoter changes to the number of observed indel events, for each pair of regulons. A promoter change occurs during the comparison of DNA sequence upstream of orthologous genes in which one gene contains a predicted binding site with a significant score while the second either does not contain a predicted binding site, or contains a predicted binding site that has a score that is not significant. The results confirmed that the predominant mechanism for divergence of these regulons relied on the assignment of new SigB promoters to genes from a common gene pool (Figure2). We extended this comparison by constructing a heat map summarizing the structures of the predicted SigB regulons in the Bc sensu lato group (Figure3; a larger version of this figure with legible annotation is presented as Additional file1: Figure S1). This heat map is color-coded: green blocks indicate the presence of a gene predicted to be in a TU controlled by a SigB promoter, while blue blocks indicate the presence of an orthologous gene not controlled by SigB. Red blocks indicate that the ortholog in question is absent from the genome.

The core SigB regulon

All of the SigB-controlled genes originally described by Van Schaik et al.[44] are found in this core regulon. We identified 27 additional genes that are predicted to be in the core SigB regulon (Figure3 and Additional file2: Table S1). Many of the additional genes encode regulatory proteins, whose altered transcription may be difficult to detect by microarray analysis, and whose regulatory control may be difficult to ascertain in the absence of experimental data. Another large group of genes added to the core regulon are predicted to encode hypothetical proteins. As a consequence, it is difficult in most instances to predict how these additional genes contribute to the function of the core SigB stress response.

The discrepancy between our predictions and the data of Van Schaik et al. is not surprising; similar difficulties in cataloguing the complete SigB-dependent response have been encountered in other systems. For instance, as noted by Nannapaneni et al.[48], three independent B. subtilis studies identified over 100 SigB-dependent target genes, but only 67 genes were commonly identified in all studies. After extensive testing in a variety of experimental conditions including ethanol and butanol shock, osmotic and oxidative stress, low-temperature growth and heat shock, Nannapaneni et al. demonstrated 166 genes that appeared to be part of the B. subtilis SigB-controlled response; 19 of these also were controlled by secondary regulators. Earlier discrepancies were likely due to differences in growth conditions, microarray platforms and the experimental setups employed. Similarly, the experiments of Van Schaik et al. likely identified a subset of the total SigB-dependent genes in Bc ATCC14579, as these experiments employed a short heat shock[44] or overexpression of the SigB protein[45], and were not as extensive as those of Nannapaneni et al. Secondary regulators interacting with the SigB regulon also are a likely complication. In S. aureus, SpoVG fine-tunes SigB-dependent regulation by antagonizing SigB[49]. In L. monocytogenes, SigB regulatory control is intertwined with that of HrcA[50], CtsR[51] and PrfA[52], complicating the analysis of SigB-mediated effects. In addition to the secondary regulators found in gene clusters 1 and 2 (Figure3), additional regulatory proteins are included in clade-specific SigB regulons (see below), and may modulate the control of gene expression by SigB, confounding the results of in vitro experiments. Thus, gene clusters 1 and 2 in Figure3 likely are a more comprehensive representation of the total SigB core regulon in the Bc sensu lato group than the early studies of Van Schaik et al..

Additions to the core SigB regulon

Although Clades C/D genomes encode cluster 5 genes, they are not functionally included in the SigB generalized stress response. Instead, Clades C/D have augmented the core SigB regulon with 24 genes from three different clusters (3a, 4a/b, in 15 TUs; Table4). Similar to gene cluster 5 in Clade A, gene clusters 3a and 4a have been added to the core SigB regulon primarily by SigB promoter assignment from a common gene pool. By contrast, gene cluster 4b was added to the SigB regulon in Clades C/D by a mixture of HGT/indels and promoter assignment (Figure3); most of these genes are not present in Clades A and B (red blocks in Figure3, cluster 4b).

Clade A and Cluster 5

Cluster 5 includes a two-component regulator not orthologous to those included in Clusters 3a/4a. This suggests that Cluster 5 adds to the core regulon a regulatory cascade that extends further into the Bc sensu lato transcriptome, and specific to Clade A. This extension of the core SigB regulon would allow the stress response to be coordinated with the environmental signal(s) to which this secondary regulator independently responds. Gene cluster 5 also appears to coordinate the enhanced synthesis of cardiolipin (CL), stimulated by at least three genes in this cluster, including CL synthase (TU #1), and phosphatidylglycerophosphate synthase and 3-ketoacyl-(acyl-carrier-protein) reductase (TU #2). Increased CL synthesis would change the composition of the cell membrane, increasing the hydrophobicity and viscosity of this membrane[54], which could have varied phenotypic effects. Increased CL levels stimulate protein translocation across the cell membrane in B. subtilis[55] and E. coli[56], and are important for high osmolarity survival in B. subtilis[57] and S. aureus[58]. High CL concentrations decrease cell envelope permeability, affecting nonspecific antibiotic resistance in E. coli[59], resistance to organic solvents[60], daptomycin resistance in S. aureus[61] and Enterococcus faecalis[62], and resistance to antimicrobial peptides such as enterococcal AS-48[63] and platelet microbicidal peptide tPMP-1[64]. CL-rich domains in B. subtilis are preferentially distributed at the medial septa and the poles during exponential growth, and the polar septal membrane and the engulfment and forespore membranes during sporulation[65]. B. subtilis spore membranes have a significantly higher CL content than membranes of exponentially growing cells, and CL appears to be essential to the proper functioning of germinant receptors[65].

Due to the extensive divergence within ABC transporter gene family, functional predictions may be difficult to substantiate by sequence similarity alone. Nevertheless, the ABC transporter genes included in TU#1 and TU#3 appear to be most similar to the ABC A superfamily[66], typically responsible for export of hydrophobic compounds. Prokaryotic members of this large gene family include the Streptomyces peucetius[67] DrrAB exporter for hydrophobic compounds such as daunomycin and doxorubicin, and the B. licheniformis BcrABC bacitracin-resistance proteins[68]. Consequently, the ABC transporters in Cluster 5 may be efflux systems either responsible for CL export or working in concert with enhanced CL production to increase resistance to deleterious hydrophobic compounds or antimicrobial peptides.

Cluster 5 TU #3 also includes the genes encoding RecA and CinA. In S. pneumoniae, the CinA ortholog appears to direct RecA to the membrane and enhance competency[69], although in B. subtilis, CinA appears to be a nucleoid-associated protein[70]. However, the most obvious consequence of enhanced RecA expression is an increased likelihood of activation of the SOS stress response by cleavage of LexA. In B. subtilis, at least 33 genes in 18 TUs appear to be directly under RecA/LexA control[71]. If the Bc sensu lato SOS response is similar, this also would extend the SigB generalized stress response into a region of the Clade A transcriptome to make these organisms more apt to express, among others, components of the excision, recombinational and error-prone repair pathways for dealing with DNA damage. Consequently, the collective inclusion of the genes in cluster 5 in the SigB-mediated stress response may fine-tune the stress response for competition with other microorganisms in a mixed microbial milieu. Cluster 5 genes also may confer some advantage during pathogenesis, for example by increasing resistance to the bactericidal effects of neutrophil-induced DNA damage or antimicrobial peptides.

Clades C/D and Clusters 3a, 4a and 4b

These clusters include the genes for three signal transduction histidine kinases, two of which are in TUs with a (probable) cognate response regulator (Table4). Two additional transcriptional regulators and a methyl-accepting chemotaxis protein also are included in these gene clusters. This suggests that an entirely different regulatory effector gene cascade extends the core SigB regulon in Clade C/D organisms, presumably coordinating the SigB generalized stress response with different environmental signals than that of Clade A organisms. Note that the two Cluster 4b regulatory proteins are unique to the genomes of invasive Clade C/D organisms, suggesting that these genes likely control functions unique to the invasive phenotype. By contrast, the SigB-responsive regulators in Clusters 3a/4a are part of the common gene pool shared with all genomes, although not included in the SigB regulon in those other organisms. This suggests that regulators found in Clusters 3a/4a are important for supervising aspects of the common metabolism shared with other Bc sensu lato organisms, but their inclusion in the generalized stress response is uniquely important for Clades C/D organisms.

Together, these gene clusters allow SigB to coordinate the expression of functions devoted to cell wall and spore structure and biosynthesis. Cluster 3a encodes an FtsI ortholog, also known as PBP-2B[72], which appears to be recruited late to the division septum in B. subtilis[73]. One of the cluster 4a regulatory proteins is the ortholog of the Bacillus LytR protein, a member of the LytR-CpsA-Psr gene family found in most members of the Firmicutes[74]. In B. subtilis, LytR lies in an operon with and controls expression of genes encoding the major autolytic amidase LytC, and a modifier gene LytB[75]. In Streptococcus pneumoniae, LytR is essential for normal septum formation[76], while the S. mutans LytR controls autolytic activity[77]. The Staphylococcus aureus LytR ortholog MsrR plays a role in cell envelope maintenance, cell separation and virulence, and appears to be connected to the sarA attenuation and virulence regulatory network[78, 79]. Thus, inclusion of LytR in the Clades C/D SigB regulons suggests that controlling cell wall-specific functions is important for supporting invasive disease.

Cluster 4b includes an S-layer protein with three SLH domains, previously designated as BslI by Kern and Schneewind[80]. While no function has yet been ascribed to BslI, the genome of Ba Sterne encodes a suite of at least 24 S-layer proteins. Eight other S-layer proteins are involved in peptidoglycan metabolism, with an additional seven important for virulence[80], suggesting that the function of the cell-surface S-layer protein BslI is similar.

Cluster 4a also encodes UDP-glucose 4-epimerase, which catalyzes the conversion of UDP-glucose to UDP-galactose[81]. Transcriptionally linked to this gene by SigB, Cluster 3a includes a member of the GNAT-transacetylase superfamily[82]. Some members of this superfamily use phosphoglucosamine as a substrate for transacetylation[82]. Collectively, these enzymes could collaborate to promote a burst in the synthesis of galactose-N-acetylglucosamine (GalNAc) during the stress response; this disaccharide is an important component of secondary cell wall polysaccharides found in Ba and some virulent Bc strains[83]. Additionally, GalNAc is an important linker molecule for tethering oligosaccharides to the collagen-like exosporium glycoprotein BclA[81], an immunodominant Ba spore antigen[84].

Cluster 3a encodes the gene for spore germination protein GerPF. In Ba spores, impairment of GerPF expression causes a germination defect in response to nutrient germinants[85]. This suggests that the SigB stress response in Clades C/D organisms has been adjusted to ensure that a sufficient amount of GerPF is made for proper spore germination potential. Targeting expression of GerPF via SigB does not appear necessary for Clades A/B organisms; possibly GerPF is more important for the proper functioning of spores produced by invasive pathogens. For instance, this protein could be important for proper germination of spores within macrophages, an event that is associated with the ability of Ba strains to cause disease[86]. Note also that some cluster 4a genes, while present in Clade D organisms, are not included in the SigB regulon (red blocks in the heat map in the cluster 4a region). This may suggest that the generalized stress response in Clade D microbes is somewhat less capable of augmenting the pathogenic potential of these organisms, in contrast to Clade C organisms, which include only virulent Ba.

Evolution of the SigB regulon and pathogenesis of the Bc sensu lato group

This sample size is admittedly small, and so the generalizations leading to these hypotheses require further exploration. Nevertheless, the hypotheses drawn from this heat map and cluster analyses are provocative. The simplest interpretation of the accumulated data (Figure4) suggests that a Clade B-like organism arose first from a B. cytotoxicus-like progenitor after receiving the SigB operon during an HGT event. Following this, a nascent core SigB regulon was assembled by the addition of SigB promoters to TUs from the B. cytotoxicus-like progenitor. Other genes were added to the core regulon by a combination of HGT and gene duplication/divergence, both events accompanied by SigB promoter assignment. Clade A and Clades C/D organisms subsequently diverged from this common Clade B-like organism by the addition of gene clusters 5 and 3a/4a/b, respectively, to the SigB core regulon. These additions occurred primarily by the assignment of SigB promoters to genes from a common gene pool, while the appearance of new genes, presumably by HGT or gene duplication, appears to be a much less common event. The TUs augmenting the SigB regulon of Clade A organisms may enhance the basal pathogenic potential of this clade, but may also make these organisms more capable of competing in the soil or other microbial consortia. Clades C/D organisms are invasive pathogens; we hypothesize that the augmentation to the SigB core regulon in these organisms is important for supporting this phenotype. Note also that some of the Clade D organisms have been referred to as ‘anthrax-like’ bacteria[27, 87], despite the tight clonality of the Ba phylogeny that excludes these organisms. The close similarity of the SigB regulon of Clade C and D organisms suggests a functional relationship between these two groups that may be important for their similar pathogenic potential.

Genome sequences used in this study

The completed genome sequences of 20 members of the B. cereus sensu lato group used in these studies are listed in Table1. In this study, we made the simplifying assumption to focus exclusively on whole genome sequences and excluded the extensive collection of draft sequences available for the Bc sensu lato group. This was due to our intent to compare the frequency of indels to promoter differences in the SigB regulons of these organisms. Lack of a particular protein coding sequence in a draft genome could be due either to deletion of this coding sequence from the genome, or a lack of sampling of this region during the sequencing for the draft assembly. A related concern was that many draft sequences likely contain regions of misassembly or low sequence quality that could complicate the analyses that we intended, since searching for specific promoter sequences in the genome would require high quality data at all points in the sequence.

Sigma factor B comparisons

To confirm the identity of the annotated SigB orthologs in each genome, all sigma factor proteins in this genome set were identified using blastp, with the Bc ATCC14579 SigB (locus tag BC1004) as the query sequence. These sequences were compared with a maximum likelihood algorithm using PhyML [Seaview[88]] with 100 bootstrap iterations and then analyzed manually to confirm the identity of all BC1004 orthologs. We further examined these proteins to confirm that the DNA binding sites encoded in each protein were the same, and were distinct from those of other sigma factors. The DNA binding sites of all sigma factor proteins encoded in the Bc ATCC14579 genome were compared using Sigma70_r2.hmm and Sigma70_r4.hmm (hidden Markov models of sigma factor region 2 and region 4, respectively), obtained from the Pfam website[89]. The utility hmmalign [hmmer-3.0b3[90]] was used to create alignments of region 2 and region 4 of all Bc ATCC14579 primary alternative sigmas, which were manually inspected. A similar comparison was made with all SigB sequences using the utility hmmalign along with Sigma70_r2.hmm and Sigma70_r2.hmm. These comparisons confirmed that all SigB orthologs encoded in these genomes had a conserved DNA binding domain, and therefore are expected to recognize a conserved DNA binding site.

SigB binding site model building and validation

We began model building by using experimentally-derived SigB binding sites derived from Bc ATCC14579[44, 45]. This data set was expanded for training purposes using phylogenetic footprinting[91]. To do this, we began with the 5’ gene of each Bc ATCC14579 operon regulated by SigB, including BC0863, BC0998, BC1001, BC1002, BC1005, BC1009, BC1010, BC3132, and BC4641. Blastp[92] was used with default parameters, except an e-value cut-off of 1e-5, to identify similar sequences in the remaining 18 Bc sensu lato genomes. We did not use the B. cytotoxicus genome in these model-building experiments, as this organism lacks a SigB operon[4]. Orthologs to the Bc ATCC14579 proteins in each organism were identified by constructing gene phylogenies for each gene family using Maximum Likelihood analysis (Seaview) and 100 bootstrap iterations. Tblastn searches were performed to ensure that any missing protein sequences in this data set were due to deletion events in the genome in question, rather than being missed during annotation. Then, the 150bp DNA sequence 5’ to each ortholog were collected and analyzed using the MEME default parameters[93] to identify the SigB binding motif 5’ to each gene. These sequences defined the initial training set used for model building. Redundancy in this training was reduced using the purge function from the MEME package. We built two complementary models of the SigB binding site to be used in parallel for all subsequent analyses. A hidden Markov model (HMM) of the training sequences was created using hmmer-1.8.5[94], while a position weight matrix (PWM) was derived using MEME. Leave-one-out cross validation[95] was used to determine an appropriate empirical p-value cut-off of .00273 for the HMM and of 4.61e-7 for the PWM. Background scores for empirical p-value calculations for the HMM analysis were determined by searching the SigB HMM against the randomly shuffled Bc ATCC14579 genome [shuffleseq from EMBOSS-6.3.1[96]]. Empirical p-values for the SigB HMM alignments were determined by comparing the alignment score to the background scores using a Perl script and R-2.13.0.

Orthologous Gene Clustering, Operon and Regulon predictions

Initially, the Ba str. Sterne transcriptome structure was extracted from the supplementary information provided by Passalacqua et al.[28]. For comparing the predicted transcriptomes from all genomes, we clustered all annotated protein-coding sequences using CD-hit[97] to identify orthologous gene clusters. This resulted in a total of 16,835 gene clusters, with cluster sizes ranging from 1 to 19 genes; each cluster was assigned a unique identifier (bID_Cluster_number; see Additional file3: Table S2). All predicted SigB-regulated genes were combined into a separate file and also clustered with CD-hit. On both data sets, we used the CD-hit default parameters except requiring an 85% sequence identity and an 80% sequence length difference cut-off. Operon predictions were performed using operonMBP[98]. Receiver operator characteristic analysis[99] was applied to the Ba str. Sterne transcriptome data to determine the appropriate cut-off to maximize true positive and minimize false positive transcription unit (TU) predictions (Figure5). The appropriate window size was determined to be 250, and the appropriate cut-off for pairing two consecutive genes as co-operonic was 14. All TU predictions were then stored in a MySQL database (see Additional file4: Table S3 for a summary of the TU predictions for each genome). To predict SigB regulon structure in each organism, the SigB HMM and PWM were used to independently scan across all 20 Bacillus genomes, and TUs 3’ to a significant binding site were extracted. We defined significant SigB binding sites as DNA regions that aligned to both the SigB HMM and PWM models, with p-values equal to or below at least one of their respective p-value cut-offs. Further, we allowed a 5’ untranslated region of up to 4.2 kb, based on the observations of Nicholas et al. in B. subtilis[100]. This information was then added to the MySQL database.

Finally, the regulon predictions were tabulated in a table of genomes vs. orthologous gene clusters and manually edited for accuracy. Scores of 0, 1, or 2 were assigned to the following states: 0, a gene not present in an orthologous cluster for a specific organism; 1, gene present but not predicted to be regulated by SigB in an orthologous cluster for a specific organism; and 2, gene present and predicted to be regulated by SigB in an orthologous cluster for a specific organism. Then, two comparisons were performed. We first used a pairwise comparison to assess the contribution of SigB promoters and indels to the structure of each regulon, and the changes between genomes. A promoter change occurs during the comparison of DNA sequence upstream of orthologous genes in which one gene contains a predicted binding site with a significant score while the second either does not contain a predicted binding site, or contains a predicted binding site that has a score that is not significant. And, we defined indels simply as the presence or lack of a protein coding sequence at a specific location; other more complex alterations, such as deletion of untranslated regions, SNPs or inversions, were not tabulated. All pairwise regulon comparisons were made in which the number of changes from 2 to 1, or 2 to 0 were counted. These numbers were normalized to the largest value counted and then plotted as indels vs. promoter changes (Figure2). Secondly, the regulon scores were placed in a table that was converted to a phylip file and analyzed using PHYLIP with Maximum Parsimony and 100 bootstrap iterations to create a SigB regulon phylogeny. This table was also analyzed in R to create a heatmap (Figure3 and Additional file1: Figure S1) using heatmap.2 from the library package gplots. Cluster analysis for this heatmap was performed using the default hclust function in R, which specifies a complete-linkage clustering method similar to neighbor-joining.