Skip to content

Advertisement

You're viewing the new version of our site. Please leave us feedback.

Learn more

BMC Evolutionary Biology

Open Access

Pan-genome and phylogeny of Bacillus cereus sensu lato

BMC Evolutionary BiologyBMC series – open, inclusive and trusted201717:176

https://doi.org/10.1186/s12862-017-1020-1

Received: 22 June 2017

Accepted: 26 July 2017

Published: 2 August 2017

Abstract

Background

Bacillus cereus sensu lato (s. l.) is an ecologically diverse bacterial group of medical and agricultural significance. In this study, I use publicly available genomes and novel bioinformatic workflows to characterize the B. cereus s. l. pan-genome and perform the largest phylogenetic and population genetic analyses of this group to date in terms of the number of genes and taxa included. With these fundamental data in hand, I identify genes associated with particular phenotypic traits (i.e., “pan-GWAS” analysis), and quantify the degree to which taxa sharing common attributes are phylogenetically clustered.

Methods

A rapid k-mer based approach (Mash) was used to create reduced representations of selected Bacillus genomes, and a fast distance-based phylogenetic analysis of this data (FastME) was performed to determine which species should be included in B. cereus s. l. The complete genomes of eight B. cereus s. l. species were annotated de novo with Prokka, and these annotations were used by Roary to produce the B. cereus s. l. pan-genome. Scoary was used to associate gene presence and absence patterns with various phenotypes. The orthologous protein sequence clusters produced by Roary were filtered and used to build HaMStR databases of gene models that were used in turn to construct phylogenetic data matrices. Phylogenetic analyses used RAxML, DendroPy, ClonalFrameML, PAUP*, and SplitsTree. Bayesian model-based population genetic analysis assigned taxa to clusters using hierBAPS. The genealogical sorting index was used to quantify the phylogenetic clustering of taxa sharing common attributes.

Results

The B. cereus s. l. pan-genome currently consists of ≈60,000 genes, ≈600 of which are “core” (common to at least 99% of taxa sampled). Pan-GWAS analysis revealed genes associated with phenotypes such as isolation source, oxygen requirement, and ability to cause diseases such as anthrax or food poisoning. Extensive phylogenetic analyses using an unprecedented amount of data produced phylogenies that were largely concordant with each other and with previous studies. Phylogenetic support as measured by bootstrap probabilities increased markedly when all suitable pan-genome data was included in phylogenetic analyses, as opposed to when only core genes were used. Bayesian population genetic analysis recommended subdividing the three major clades of B. cereus s. l. into nine clusters. Taxa sharing common traits and species designations exhibited varying degrees of phylogenetic clustering.

Conclusions

All phylogenetic analyses recapitulated two previously used classification systems, and taxa were consistently assigned to the same major clade and group. By including accessory genes from the pan-genome in the phylogenetic analyses, I produced an exceptionally well-supported phylogeny of 114 complete B. cereus s. l. genomes. The best-performing methods were used to produce a phylogeny of all 498 publicly available B. cereus s. l. genomes, which was in turn used to compare three different classification systems and to test the monophyly status of various B. cereus s. l. species. The majority of the methodology used in this study is generic and could be leveraged to produce pan-genome estimates and similarly robust phylogenetic hypotheses for other bacterial groups.

Keywords

Bacillus cereus sensu lato Bacillus cereus group Bacillus Pan-genomePhylogenyPhylogeneticsPan-GWASBayesian model-based clusteringPhylogenetic clustering

Background

Bacillus cereus sensu lato (s. l.) is an ecologically diverse bacterial group that comprises a growing number of species, many of which are medically or agriculturally important. Historically recognized and most well-sampled of the species are B. anthracis (the causative agent of anthrax), B. cereus sensu stricto (capable of causing food poisoning and other ailments), and B. thuringiensis (used to control insect pests). Other species are distinguished by rhizoidal growth patterns (B. mycoides and B. pseudomycoides [1]), thermotolerance and cytotoxicity (B. cytotoxicus [2]), psychrotolerance and ability to cause food spoilage (B. weihenstephanensis [3] and B. wiedmannii [4]), and utility as a probiotic in animal nutrition (B. toyonensis [5]). In addition, several new species have also recently been described (B. bingmayongensis [6], B. gaemokensis [7], and B. manliponensis [8]). In order to understand the fantastic diversity of B. cereus s. l. and its concomitant ability to occupy diverse environmental niches and exhibit a variety of phenotypes, it is crucial to accurately characterize genomic diversity within the group and to generate robust phylogenetic hypotheses about the evolutionary relationships among group members.

A typical B. cereus s. l. genome contains ≈5500 protein-coding genes [9, 10]. Due to rampant horizontal gene transfer in bacterial ecosystems, however, the genome of a particular strain or species often contains genes not found in closely related taxa [11]. Thus, it is now common practice to characterize the full gene complement of a closely related group of bacterial strains or species, otherwise known as a “pan-genome” [11]. In this study, a “core” gene is defined as present in at least 99% of sampled taxa, an “accessory” gene as a non-core gene present in at least two taxa, and a “unique” gene as present in only one taxon. A previous effort to characterize the B. cereus s. l. pan-genome [12] based on a comparison of a relatively small number of strains estimated that there are ≈3000 core genes and ≈22,500 total genes in the B. cereus s. l. pan-genome. A more recent study [13] using 58 strains reported similar estimates.

Phylogenetic hypotheses of B. cereus s. l. have been generated from a variety of data sources, including 16S rRNA sequences [12], amplified fragment length polymorphism (AFLP) data [14, 15], multilocus sequence typing (MLST) of housekeeping genes [1518], single-copy protein-coding genes [19], locally collinear blocks (LCBs) [13], conserved protein-coding genes [13], whole-genome single nucleotide polymorphisms (SNPs) [18], and digital DNA-DNA hybridization (dDDH) data [20]. Phylogenetic analyses have used distance methods [13, 14, 17, 20], maximum likelihood [14, 18, 19], maximum parsimony [14], and Bayesian methods [16]. For the most part, published phylogenies have tended to agree with and reinforce one another, although naturally there have been different classification systems developed with attendant implications for species designations. One popular classification system divides the B. cereus s. l. phylogeny into three broad clades [13, 16, 21]; traditionally, Clade 1 contains B. anthracis, B. cereus, and B. thuringiensis; Clade 2 contains B. cereus and B. thuringiensis; and Clade 3 contains a greater diversity of species including B. cereus, B. cytotoxicus, B. mycoides, B. thuringiensis, B. toyonensis, and B. weihenstephanensis. A somewhat more fine-grained classification system divides the phylogeny into seven major groups [14, 15, 18], each with its own thermotolerance profile [14] and propensity to cause food poisoning [22].

In this study I aimed to produce the most accurate and comprehensive estimate of the B. cereus s. l. pan-genome and phylogeny to date by analyzing all publicly available B. cereus s. l. genome data with a novel bioinformatic workflow for pan-genome characterization and pan-genome-based phylogenetic analysis.

Methods

Distance-based phylogeny of the genus Bacillus

All “reference” and “representative” Bacillus genome assemblies were retrieved from the NCBI RefSeq [23] database in October 2016, comprising 86 assemblies from 74 well-described Bacillus species and 44 assemblies from as-yet uncharacterized species. In addition, 16 “latest” assemblies were added for five Bacillus species that are thought to be part of B. cereus s. l. (B. bingmayongensis [6], B. gaemokensis [7], B. pseudomycoides [1], B. toyonensis [5], and B. wiedmannii [4]). In total, 146 Bacillus genomes were included in the distance-based phylogenetic analysis. The sketch function in Mash [24] version 1.1.1 (arguments: -k 21 -s 1000) was used to create a compressed representation of each genome, and then the Mash distance function was used to generate all pairwise distances among genomes. The Mash distance matrix was converted to PHYLIP format and analyzed with FastME [25] version 2.1.4 using the default BIONJ [26] algorithm.

Creation of taxon sets

BCSL_114

All complete genomes of eight B. cereus s. l. species (B. anthracis, B. cereus, B. cytotoxicus [2], B. mycoides, B. pseudomycoides [1], B. thuringiensis, B. toyonensis [5], and B. weihenstephanensis [3]) were downloaded from the NCBI RefSeq [23] database in October 2016, which altogether comprised 114 genomes. One strain from each species was designated the “reference taxon” for that species, as required by HaMStR [27] (Table 1). This taxon set of complete genomes (“BCSL_114”; Table 1) was used to build the HaMStR databases and as the basis for the majority of the analyses performed in this study.
Table 1

Species composition of taxon sets

Species

Representatives in bcsl_114

Representatives in bcsl_498

HaMStR reference taxon

RefSeq assembly accession

B. anthracis

42

128

Ames

GCF_000007845

B. bingmayongensis

0

1

  

B. cereus (sensu stricto)

30

258

ATCC 14579

GCF_000007825

B. cytotoxicus

1

2

NVH 391-98

GCF_000017425

B. gaemokensis

0

2

  

B. manliponensis

0

1

  

B. mycoides

2

13

ATCC 6462

GCF_000832605

B. pseudomycoides

1

1

DSM 12442

GCF_000161455

B. thuringiensis

35

73

97-27

GCF_000008505

B. toyonensis

1

1

BCT-7112

GCF_000496285

B. weihenstephanensis

2

6

KBAB4

GCF_000018825

B. wiedmannii

0

11

  

Bacillus sp.

0

1

  

BCSL_498

To perform analyses involving all publicly available B. cereus s. l. genome data, all “latest assemblies” were downloaded for the eight species mentioned above, and based on analysis of the Bacillus distance-based phylogeny, assemblies were added for B. bingmayongensis [6], B. gaemokensis [7], B. manliponensis [8], B. wiedmannii [4], and one uncharacterized species (Bacillus sp. UNC437CL72CviS29), which altogether comprised 498 genomes (“BCSL_498”; Table 1). A list of RefSeq assembly accessions for all taxa used in this study is provided (Additional file 1).

Isolate metadata

B. cereus s. l. isolate metadata, including “Assembly Accession”, “Disease”, “Host Name”, “Isolation Source”, “Motility”, and “Oxygen Requirement” was downloaded from PATRIC [28] in December 2016. This metadata was used to associate patterns of gene presence and absence with phenotypes exhibited by groups of taxa.

Genome annotation

All B. cereus s. l. genomes were annotated de novo with Prokka [29] version 1.12-beta (arguments: –kingdom Bacteria –genus Bacillus).

Pan-genome inference

The pan-genome of B. cereus s. l. was inferred with Roary [30] version 3.7.0. The BCSL_114 Prokka annotations were provided to Roary as input; in turn, Roary produced a gene presence/absence matrix (Additional file 2), a multi-FASTA alignment of core genes using PRANK [31] version 0.140603, and a tree based on the presence and absence of accessory genes among taxa using FastTree 2 [32] version 2.1.9. The “accessory binary tree” was computed using only the first 4000 genes in the accessory genome.

Phylogenetic network analysis

A NEXUS-format binary version of the BCSL_114 gene presence/absence matrix was analyzed with SplitsTree4 [33] version 4.14.4. Three methods of calculating distances between taxa were evaluated: Uncorrected_P, GeneContentDistance [34], and the MLDistance variant of GeneContentDistance [34]. The NeighborNet [35] algorithm was used to reconstruct the phylogenetic network.

Genotype-phenotype association

Scoary [36] version 1.6.9 was used to associate patterns of gene presence and absence with particular phenotypes (traits), an analysis known as “pan-GWAS” [36]. Scoary required two basic input files: the BCSL_114 gene presence/absence matrix, augmented with gene presence/absence information for BCSL_498 taxa obtained from orthology determination with HaMStR [27] (Additional file 3), and a binary trait matrix that was created using the isolate metadata obtained from PATRIC (Additional file 4). Assignment of traits to taxa was performed conservatively in that missing data was not assumed to be an indication of the presence or absence of a particular trait. Scoary was run with 1,000 permutation replicates, and genes were reported as significantly associated with a trait if they attained a naive P-value less than 0.05, a Benjamini-Hochberg-corrected P-value less than 0.05, an empirical P-value less than 0.05, and were not annotated as “hypothetical proteins”. Lists of genes were subsequently tested for enrichment of biological processes using the data and services provided by AmiGO 2 [37] version 2.4.24, which in turn used the PANTHER database [38] version 11.1.

HaMStR database creation

The orthologous protein sequence clusters output by Roary were filtered to produce a set of gene models suitable for use with HaMStR [27] version 13.2.6. HaMStR enables one to build gene models for a clade of interest (using, ideally, high-quality complete genomes), which are subsequently used to identify orthologs in other sequence data (e.g., draft genome assemblies, transcriptomes, etc.). HaMStR required that each sequence cluster contain at least one sequence from the set of previously selected reference taxa (Table 1), so clusters not meeting this requirement were omitted. Furthermore, each cluster was required to contain at least four sequences, the minimum number of sequences required to produce an informative unrooted phylogenetic tree. Finally, clusters that Roary flagged as having a quality-control issue were removed. The 9,070 clusters that passed these filters were aligned using the linsi algorithm in MAFFT [39] version 7.305. Gene models (i.e., Hidden Markov Models, or HMMs) were produced from the aligned cluster sequences using the hmmbuild program from HMMER [40] version 3.0. Finally, for each reference taxon, a BLAST [41] database was built using the full complement of protein-coding genes for that taxon. This completed the construction of the initial HaMStR database, which is called “HAMSTR_FULL”. A variant of HAMSTR_FULL called “HAMSTR_CORE” was created, which contained only the 594 gene models corresponding to core genes.

Mobile genetic element removal

For tree-based phylogenetic analyses that assume a process of vertical inheritance, the inclusion of mobile genetic elements (MGEs) that may be horizontally transferred is likely to confound the phylogenetic inference process [42]. Thus, an effort was made to identify and remove putative MGEs from the HaMStR databases. In December 2016, a list of Bacillus genes derived from a plasmid source was downloaded from the NCBI Gene [43] database. In addition, all genes were exported from the ACLAME [44] database version 0.4. Using this information, gene models that were either plasmid-associated or found in the ACLAME list of MGEs were removed from HaMStR databases. Gene models whose annotation included the keywords “transposon”, “transposition”, “transposase”, “insertion”, “insertase”, “plasmid”, “prophage”, “intron”, “integrase”, or “conjugal” were also removed. The resulting HaMStR databases, “HAMSTR_FULL_MGES_REMOVED” and “HAMSTR_CORE_MGES_REMOVED”, contained 8954 and 578 gene models, respectively. The workflow used to construct the HAMSTR_FULL_MGES_REMOVED database is shown as a diagram (Additional file 5).

Orthology determination

The protein-coding gene annotations of “query” taxa — i.e., taxa not included in BCSL_114 — were searched for sequences matching HaMStR database gene models using HaMStR [27] version 13.2.6 (which in turn used GeneWise [45] version 2.4.1, HMMER [40] version 3.0, and BLASTP [41] version 2.2.25+). In the first step of the HaMStR search procedure, the hmmsearch program from HMMER was used to identify translated substrings of protein-coding sequence that matched a gene model in the database, which were then provisionally assigned to the corresponding sequence cluster. To reduce the number of highly divergent, potentially paralogous sequences returned by this initial search, the E-value cutoff for a “hit” was set to 1e-05 (the HaMStR default was 1.0). In the second HaMStR step, BLASTP was used to compare the hits from the HMM search against the proteome of the reference taxon associated with that gene model; sequences were only retained if the reference taxon protein used in the construction of the gene model was also the best BLAST hit. The E-value cutoff for the BLAST search was set to 1e-05 (the HaMStR default was 10.0).

Data matrix construction

Amino acid sequences assigned to orthologous sequence clusters were aligned using MAFFT [39] version 7.305. The resulting amino acid alignments were converted to corresponding nucleotide alignments using a custom Perl script that substituted for each amino acid the proper codon from the original coding sequence. Initial orthology assignment may sometimes result in multiple sequences for a particular taxon/locus combination [46], which need to be reduced to a single sequence for inclusion in phylogenetic data matrices. For this task the “consensus” [47] procedure was used, which collapsed all sequence variants into a single sequence by replacing multi-state positions with nucleotide ambiguity codes. Following application of the consensus procedure, individual sequence cluster alignments were concatenated, adding gaps for missing data as necessary using a custom Perl script. The workflow used for orthology determination and data matrix construction is shown as a diagram (Additional file 6).

Maximum likelihood phylogenetic analysis

Concatenated nucleotide data matrices were analyzed under the maximum likelihood criterion using RAxML [48] version 8.2.8 (arguments: -f d -m GTRGAMMAI). The data were analyzed either with all nucleotides included in a single data subset (ALL_NUC), or with sites partitioned by codon position (CODON_POS). Partitioned analyses assigned a unique instance of the substitution model to each data subset, with joint branch length optimization. Analyses of the BCSL_114 taxon set consisted of an adaptive best tree search [49] and an adaptive bootstopping procedure that used the autoMRE RAxML bootstopping criterion [50]; thus, the number of search replicates performed varied from 10 to 1000, depending on the analysis. DendroPy [51] was used to map bootstrap probabilities onto the best tree. Analysis of the BCSL_498 taxon set required ≈256 GB of RAM and multiple weeks of runtime, and was thus limited to a single best tree search.

Recombination detection

Genomic regions that may have been involved in past recombination events should be excluded from phylogenetic analyses that assume a process of vertical inheritance, or phylogenetic inference methods should incorporate this information to produce a more accurate phylogeny [42]. In this study, two different software packages that address this problem were evaluated. First, the profile program from PhiPack [52] was used to flag and remove from concatenated data matrices sites that exhibited signs of mosaicism. Following the procedure employed in Parsnp [53], the profile program defaults were used, except that the step size was increased from 25 to 100 (-m 100). RAxML was then used to create new versions of data matrices that excluded regions whose Phi statistic P-value was less than 0.01. Second, ClonalFrameML [54] (downloaded from GitHub June 14, 2016) with default parameters was used to correct the branch lengths of phylogenies to account for recombination. ClonalFrameML required all ambiguous bases in data matrices to be coded as “N”.

Maximum parsimony phylogenetic analysis

Concatenated nucleotide data matrices were analyzed under the maximum parsimony criterion using PAUP* [55] version 4.0a150. A heuristic search was performed using default parameters.

Tree distance calculation

To quantify the difference between pairs of tree topologies, both the standard and normalized Robinson-Foulds distance [56] were calculated with RAxML [48] version 8.2.8 (arguments: -f r -z).

Tree visualization

All visualizations of phylogenetic trees were produced with FigTree [57] version 1.4.2, except Fig. 4, which was produced with iTOL [58] version 3.5.3.

Taxon clustering

To complement phylogenetic analysis and existing classification systems, taxa were clustered with hierBAPS [59] (bugfixed version dated August 15, 2013), a Bayesian model-based population genetic approach that accounts for admixture within and among lineages. The BCSL_498 alignment of 8954 genes (MAT_6) was provided to hierBAPS as input, and hierBAPS was directed to produce a single-level clustering with a maximum of 10 clusters.

Clustering of taxon-associated attributes

The degree of clustering of taxa sharing a common attribute, given a phylogeny relating those taxa, was quantified using the genealogical sorting index [60] (gsi) version 0.92 made available through the web service at molecularevolution.org [61]. Significance of the gsi was determined by running 104 permutation replicates.

Results

Distance-based phylogeny of the genus Bacillus

The Mash-distance-based phylogeny of the genus Bacillus (Additional file 7) indicated a B. cereus s. l. clade containing the following species: B. anthracis, B. bingmayongensis, B. cereus (sensu stricto), B. cytotoxicus, B. gaemokensis, B. manliponensis, B. mycoides, B. pseudomycoides, B. thuringiensis, B. toyonensis, B. weihenstephanensis, B. wiedmannii, and one uncharacterized species (Bacillus sp. UNC437CL72CviS29). Within B. cereus s. l., the first taxon to split off from the remainder of the group was B. manliponensis, followed by B. cytotoxicus (which has been previously recognized as an outlier [12, 19]).

Pan-genome inference

The pan-genome of B. cereus s. l. was inferred with Roary [30] using the BCSL_114 taxon set. Roary produced a total of 59,989 protein-coding gene sequence clusters (Additional file 2) from an average of 5726 genes per input genome (Additional file 8). The shortest cluster sequence was 122 nt, the longest cluster sequence was 22,967 nt, and the average length of a cluster sequence was 755 nt. The average difference between the shortest and longest sequence in a cluster was only 67 nt, suggesting that the input data was relatively uniform, as would be expected from complete genomes. The B. cereus s. l. “core genome”, consisting of genes present in at least 99% of taxa sampled, was represented by 598 genes (≈1% of all genes). A rarefaction curve shows that after ≈35 genomes have been sampled (≈31% of all genomes), the number of core genes remains fairly constant at ≈600 genes, while the total number of genes in the pan-genome continues to increase almost linearly (Additional file 9). The 59,391 non-core genes were divided into 32,324 “accessory genes” (i.e., non-core genes present in at least two taxa; ≈54% of all genes), and 27,067 “unique genes” (i.e., genes present in only one taxon; ≈45% of all genes). A rarefaction curve shows that as genomes are sampled, genes never before observed continue to be found at a fairly steady rate, and the total number of unique genes discovered continues to increase, with no indication of soon approaching an asymptote (Additional file 10); these trends indicate that the B. cereus s. l. pan-genome is “open”. Finally, Roary produced an “accessory binary tree”, which was plotted alongside core and accessory gene presence/absence information (Additional file 11). This figure shows that the outermost B. cereus s. l. clades include taxa with relatively few accessory genes included in the analysis (≤40), such as B. cytotoxicus, B. mycoides, and B. pseudomycoides; by contrast, the genomes with the most accessory genes present (>1000) belong to the highly clonal clade of B. anthracis strains.

Phylogenetic network analysis

SplitsTree4 [33] was used to build a phylogenetic network from the gene presence/absence information produced by Roary. The choice of method for computing distance did affect network branch lengths; the network presented here was computed with the MLDistance variant of GeneContentDistance (Fig. 1), as that method seemed most appropriate for gene presence/absence data. The phylogenetic network recapitulated both the three-clade [13, 16, 21] and seven-group [14, 15, 18] classification systems used in previous studies. Group I and Group VII, both part of Clade 3, were most radically diverged from the remainder of the network. Notably, Group II was absent from the network, as it was not represented by any complete B. cereus s. l. genomes at the time the study was performed.
Fig. 1

Phylogenetic network analysis of BCSL_114. Gene presence/absence information produced by Roary was provided as input to SplitsTree, which used the MLDistance variant of GeneContentDistance together with the NeighborNet algorithm to reconstruct the phylogenetic network. Major B. cereus s. l. clades and groups are indicated, along with representative taxa. Please see Additional file 17 for high resolution image

Genotype-phenotype association

Scoary [36] was used to associate patterns of gene presence and absence with particular phenotypes (traits), an analysis known as “pan-GWAS” [36]. Pan-GWAS was performed for the following traits: isolation source (cattle, human, invertebrate, non-primate mammal, or soil); motility; oxygen requirement (aerobic or facultative); and disease (anthrax or food poisoning). Eight of ten traits tested had some number of significant positively or negatively associated genes (Table 2). Traits with a sufficient number of associated genes were tested for possible enrichment of gene ontology biological processes (Additional file 12). The most interesting findings from this analysis concerned taxa isolated from soil. Specifically, metabolic and biosynthetic processes involving quinone (and in particular, menaquinone) were positively associated with soil isolates. Analysis of quinone species present in soil have been used previously to characterize soil microbiota [62]. Furthermore, a high ratio of menaquinone to ubiquinone (the two dominant forms of quinone in soil) has been associated with the presence of gram-positive bacteria such as Bacillus species [63]. On the other hand, biological processes involving flagella, cilia, or motility more generally were negatively associated with soil isolates. This finding is consistent with observations that motility may not be necessary for bacterial colonization of plant roots [64], doubts about the evolutionary advantage of maintaining flagella in a soil environment [65], and general properties of soil that bring into question the importance of active movement and the extent to which it occurs [66].
Table 2

Scoary result summary

Trait

Taxa

Taxa

Significant

Significant

 

with

without

positively

negatively

 

trait

trait

associated

associated

   

genes

genes

Isolation source: cattle

25

331

358

227

Isolation source: human

44

312

53

47

Isolation source: invertebrate

16

340

0

0

Isolation source: non-primate mammal

46

310

162

85

Isolation source: soil

121

235

34

34

Motility

88

18

0

0

Oxygen requirement: aerobic

65

41

15

18

Oxygen requirement: facultative

41

65

20

11

Disease: anthrax

85

63

44

1

Disease: food poisoning

43

105

3

23

Concatenated data matrices

In total, six different concatenated nucleotide data matrices were constructed and analyzed (MAT_1MAT_6; Table 3). The majority of the data matrices used the BCSL_114 taxon set (MAT_1MAT_5); only MAT_6 used the BCSL_498 taxon set. Various gene sets were used, including
  1. 1)
    all core genes identified by Roary (all_core);
    Table 3

    Concatenated data matrix statistics

    Matrix

    Taxon set

    Gene set

    Alignment method

    Alignment length (nt)

    Matrix completeness

    Ambiguous characters

    mat_1

    bcsl_114

    all_core

    prank

    508,158

    97.9%

    0.4%

    mat_2

    bcsl_114

    hamstr_core

    mafft

    502,005

    98.4%

    0.0%

    mat_3

    bcsl_114

    hamstr_core_mges_removed

    mafft

    486,546

    98.4%

    0.0%

    mat_4

    bcsl_114

    hamstr_core_mges_removed + PhiPack sites removed

    mafft

    134,225

    98.0%

    0.0%

    mat_5

    bcsl_114

    hamstr_full_mges_removed

    mafft

    7,962,138

    47.4%

    0.0%

    mat_6

    bcsl_498

    hamstr_full_mges_removed

    mafft

    8,207,628

    66.1%

    0.1%

     
  2. 2)

    only the core genes used to build the HaMStR database (hamstr_core);

     
  3. 3)

    HaMStR core genes with mobile genetic elements (MGEs) removed (hamstr_core_mges_removed), and a variant of this gene set with PhiPack sites removed; and finally,

     
  4. 4)

    all HaMStR genes with MGEs removed (hamstr_full_mges_removed).

     

Aligned data matrices ranged from 96,802 nt to 8,207,628 nt in length. Matrix completeness, defined as the percentage of non-missing data, ranged from 47.4 to 99.5%. The percentage of ambiguous characters present in data matrices ranged from 0.0 to 17.0%.

Phylogenetic analyses

In total, nine different phylogenetic analyses of the six concatenated data matrices were performed (Table 4). Eight of the nine analyses used maximum likelihood (ML_1ML_8), and one analysis used maximum parsimony (MP_1). For reasons of computational tractability, all exploratory analyses used the BCSL_114 taxon set (ML_1ML_7 and MP_1); only when the best-performing methods were established was analysis of the BCSL_498 taxon set pursued (ML_8). During the exploratory phase, several variables were tested for their effect on phylogenetic outcome:
  1. 1)
    use of MAFFT instead of PRANK to align protein sequence clusters;
    Table 4

    Phylogenetic analysis statistics 1

    Analysis

    Matrix

    Partitioning scheme

    Unique alignment

    Best tree

    Bootstrap

    Nodes with

    Nodes with

    Nodes with

       

    patterns

    searches

    replicates

    B P=1.0

    B P≥0.5

    B P<0.5

    ml_1

    mat_1

    all_nuc

    46,395

    1000

    100

    68

    98

    14

    ml_2

    mat_2

    all_nuc

    46,174

    1000

    100

    68

    93

    19

    ml_3

    mat_3

    all_nuc

    44,750

    1000

    200

    63

    94

    18

    ml_4

    mat_3

    codon_pos

    49,889

    1000

    200

    64

    92

    20

    ml_5

    mat_4

    all_nuc

    11,729

    1000

    200

    53

    86

    26

    ml_6

    mat_5

    all_nuc

    691,147

    10

    100

    92

    112

    0

    ml_7

    mat_5

    codon_pos

    852,707

    28

    100

    89

    112

    0

    ml_8

    mat_6

    all_nuc

    3,948,459

    1

    0

    n/a

    n/a

    n/a

    mp_1

    mat_3

    all_nuc

    83,3832

    n/a

    n/a

    n/a

    n/a

    n/a

    1 ml = maximum likelihood; mp = maximum parsimony; BP = bootstrap probability; n/a = not applicable

    2Number of parsimony-informative characters

     
  2. 2)

    removal of MGEs;

     
  3. 3)

    use of maximum parsimony in addition to maximum likelihood;

     
  4. 4)

    partitioning of sites by codon position;

     
  5. 5)

    removal of sites implicated in recombination; and finally,

     
  6. 6)

    use of all eligible genes from the pan-genome versus only core genes.

     

Importantly, all phylogenetic analysis results recapitulated the three-clade [13, 16, 21] and seven-group [14, 15, 18] classification systems of previous studies. Taxa were consistently assigned to the same clade and group, independent of the particular phylogenetic analysis performed. Thus, topological differences between analysis results, as measured by the Robinson-Foulds distance [56] (Additional file 13), were confined to intra-group relationships. Bootstrap support was fairly consistent for all analyses that used core genes, and increased dramatically when all eligible genes from the pan-genome were used (Table 4). Additional detail about the phylogenetic analyses, and the logic behind their progression, is provided in the subsections that follow.

Choice of multiple sequence alignment program

Roary produced multiple sequence alignments of all 598 core genes with PRANK [31], which explicitly models insertions and deletions, but as a consequence runs more slowly than some other alignment programs. The PRANK alignments were concatenated to produce MAT_1. A similar matrix was built using the 594 HaMStR-eligible core genes, except that the gene sequence clusters were aligned with MAFFT [39] (MAT_2). Phylogenetic analyses of these two matrices with RAxML [48] revealed only negligible differences in bootstrap probabilities (ML_1 vs. ML_2; Table 4), so for the sake of computational efficiency MAFFT was used for the remainder of the analyses.

Removal of mobile genetic elements

For tree-based phylogenetic analyses that assume a process of vertical inheritance, the inclusion of mobile genetic elements (MGEs) that may be horizontally transferred is likely to confound the phylogenetic inference process [42]. Thus, putative MGEs were identified and removed from HAMSTR_CORE, leaving 578 core genes (HAMSTR_CORE_MGES_REMOVED). Phylogenetic analysis of this slightly smaller data matrix (MAT_3) revealed comparable bootstrap probabilities to those from the analysis that used HAMSTR_CORE (ML_3 vs. ML_2; Table 4); nevertheless, out of principle, HaMStR databases with MGEs removed were used for the remainder of the analyses.

Partitioning of sites by codon position

It is well known that nucleotides in different codon positions (first, second, or third) are likely to be under different selective pressures [67]; thus, when analyzing protein-coding nucleotide sequences, it is common practice to apply a different substitution model (or different instance of the same substitution model) to the sites associated with each codon position, thus effectively partitioning the data matrix into three data subsets. The effect of partitioning by codon position was tested with two different matrices (MAT_3 and MAT_5); only negligible differences in bootstrap probabilities were found as compared to the unpartitioned results (ML_4 vs. ML_3 and ML_7 vs. ML_6; Table 4).

Removal of sites implicated in recombination

Genomic regions that may have been involved in past recombination events should not be used for phylogenetic analyses that assume a process of vertical inheritance [42]. The profile program from PhiPack [52] was used to flag and remove sites from MAT_3 that exhibited signs of mosaicism. The resulting data matrix (MAT_4) contained less than one-fourth the number of unique alignment patterns of MAT_3, thus representing a substantial reduction in data suitable for phylogenetic analysis. This was reflected in bootstrap probabilities, which were somewhat depressed overall (ML_5 vs. ML_3; Table 4). It was thus concluded that removing sites implicated in recombination had a deleterious effect on phylogenetic analysis results, and so this procedure was not applied to subsequent analyses.

Use of all eligible genes from the pan-genome versus only core genes

Using all eligible genes (HAMSTR_FULL_MGES_REMOVED) for phylogenetic analysis as opposed to using only core genes (HAMSTR_CORE_MGES_REMOVED) caused bootstrap probabilities to increase dramatically (ML_6 vs. ML_3 and ML_7 vs. ML_4; Table 4). Thus, the ML_6 result was selected as the best estimate of the phylogenetic relationships among the BCSL_114 taxa. ClonalFrameML [54] was used to correct the branch lengths of this tree to account for recombination, and the tree was rooted using B. cytotoxicus [12, 19]. The resulting BCSL_114 phylogeny is shown as a phylogram with major clades and groups indicated (Fig. 2), and as a cladogram with bootstrap probabilities annotated (Additional file 14).
Fig. 2

BCSL_114 maximum likelihood phylogenetic analysis results. Phylogram depicting the best estimate of the phylogenetic relationships among BCSL_114 taxa, computed with RAxML using 8954 genes (ML_6; Table 4). ClonalFrameML was used to correct the branch lengths of the tree to account for recombination, and B. cytotoxicus was used to root the tree. Major B. cereus s. l. clades and groups are indicated. Please see Additional file 18 for high resolution image

Maximum likelihood-based analysis of all taxa

Once the exploratory analyses were completed, an analysis of BCSL_498 was executed using the HAMSTR_FULL_MGES_REMOVED gene set. The average number of genes included in the analysis for each species, clade, and group is given in Additional file 8, and a count of the number of genes included for each taxon is given in Additional file 15. Due to the size of the data matrix (almost 4×106 unique alignment patterns), only a single best tree search replicate was completed (ML_8; Table 4). Informed by the distance-based analysis of Bacillus species (Additional file 7), the tree was rooted using B. manliponensis. The resulting BCSL_498 phylogeny is shown as a phylogram with major clades and groups indicated (Fig. 3). In contrast to analyses of BCSL_114, Group II is now represented, and is located on the tree where expected [14, 15, 18]. Based on this topology of currently sequenced genomes, a count of the number of taxa by species is provided for major B. cereus s. l. clades and groups (Additional file 8).
Fig. 3

BCSL_498 maximum likelihood phylogenetic analysis results. Phylogram depicting an estimate of the phylogenetic relationships among BCSL_498 taxa, computed with RAxML using 8954 genes (ML_8; Table 4). B. manliponensis was used to root the tree. Major B. cereus s. l. clades and groups are indicated. Please see Additional file 19 for high resolution image

Taxon clustering

The hierBAPS [59] clustering analysis divided BCSL_498 into nine clusters (Additional file 15), which are displayed alongside major B. cereus s. l. clades and groups in Fig. 4. The hierBAPS clusters are congruent with the three-clade classification system, and largely agree with the seven-group classification system, with the following differences. Clade 1 included members of three clusters (as opposed to only two groups), and Clade 3 included members of six clusters (as opposed to four groups). Some Clade 3 clusters expanded slightly relative to their counterpart group to include taxa that were not assigned to any group, and B. manliponensis was assigned to its own cluster. Interestingly, two Clade 3 B. cereus taxa were assigned to Cluster 3, whereas other members of Cluster 3 were assigned to Clade 1, thus suggesting some genetic admixture between these two clades that was not reflected in the phylogenetic analysis.
Fig. 4

Phylogeny showing assignment of taxa to clades, groups, and clusters. Circular phylogram depicting an estimate of the phylogenetic relationships among BCSL_498 taxa, computed with RAxML using 8954 genes (ML_8; Table 4). B. manliponensis was used to root the tree. Major B. cereus s. l. clades are indicated by highlighting of taxon labels, major groups are indicated by the inner colored strip, and hierBAPS clusters are indicated by the outer colored strip. Please see Additional file 20 for high resolution image

Clustering of taxon-associated traits

The genealogical sorting index [60] (gsi) was used to quantify the degree of clustering of taxa sharing a common attribute given a phylogeny relating those taxa. The gsi statistic for a particular attribute takes a value from the unit interval [0,1]; if taxa associated with the attribute form a monophyletic group, the g s i=1; otherwise, the greater the degree to which taxa associated with the attribute are dispersed throughout the tree (accounting for the number of taxa and the size of the tree), the smaller the gsi will be for that attribute.

Quantifying the degree of B. cereus s. l. species monophyly

The gsi was calculated for six B. cereus s. l. species that were sufficiently represented in the BCSL_498 phylogeny; all P-values were < <0.05 (Additional file 16 and Table 5). Due to its highly clonal nature, B. anthracis was the species closest to monophyly (g s i=0.95), and would have indeed been monophyletic except that one B. anthracis taxon (GCF_001029875) did not group with the others (but still placed in Group III). This might be a misannotation and should be investigated. B. weihenstephanensis was the species furthest from monophyly (g s i=0.15), primarily because it was represented by only six taxa, one of which (GCF_000518025) was found in Group IV — the remainder were found in Group VI. Again, the annotation of the Group IV taxon with regard to species affiliation should be scrutinized.
Table 5

Monophyly status of B. cereus s. l. species, as quantified by the gsi

Species

Representatives in bcsl_498

gsi value

P-value

B. anthracis

128

0.95

0.0001

B. cereus

258

0.55

0.0001

B. mycoides

13

0.34

0.0001

B. thuringiensis

73

0.36

0.0001

B. weihenstephanensis

6

0.15

0.0021

B. wiedmannii

11

0.58

0.0001

Quantifying the degree of clustering of taxa sharing common traits

The gsi was calculated for ten traits shared by various B. cereus s. l. taxa using the BCSL_114 phylogeny from the ML_6 analysis; all P-values with the exception of one were less than 0.05 (Table 6). As not all of the taxa in BCSL_114 were assayed for each trait, the gsi values were artificially depressed; nevertheless, their relative values may be compared. The traits with the largest gsi values were “isolation source: cattle” and “isolation source: non-primate mammal”, the taxa associated with the former being a subset of the taxa associated with the latter. These taxa were all located in Group III, and all but two were identified as B. anthracis. This finding is consistent with the prevalence of mortality due to anthrax among cattle and other herbivores [68].
Table 6

Degree of clustering of taxa sharing common traits, as quantified by the gsi

Trait

Taxa

gsi

P-value

 

with trait

value

 

Isolation source: cattle

9

0.33

0.0001

Isolation source: human

17

0.19

0.0273

Isolation source: invertebrate

5

0.16

0.0462

Isolation source: non-primate mammal

15

0.37

0.0001

Isolation source: soil

17

0.27

0.0002

Motility

25

0.16

0.2100

Oxygen requirement: aerobic

15

0.19

0.0268

Oxygen requirement: facultative

9

0.22

0.0037

Disease: anthrax

11

0.26

0.0008

Disease: food poisoning

8

0.17

0.0302

Discussion

I show that the B. cereus s. l. pan-genome is still “open” (Additional files 9 and 10), thus implying that continued sampling of the group — especially of underrepresented taxa such as environmental strains [17] — will continue to reveal novel gene content. My estimate of the number of protein-coding genes in the B. cereus s. l. core and pan-genome (≈600 and ≈60,000, respectively), based on 114 complete genomes, is consistent with previous estimates [12, 13], as more extensive sampling of an open pan-genome will necessarily reduce the core genome size while simultaneously increasing the pan-genome size. It is interesting to observe that the basic phylogenetic structure of B. cereus s. l. can be accurately computed by relatively quick phylogenetic analyses based solely on the distribution of accessory genes among taxa (Fig. 1 and Additional file 11), which may in fact be sufficient for some applications. The diversity and adaptability of B. cereus s. l. may be in part attributable to the significant proportion of unique genes in its pan-genome (≈27,000, almost 50% of all genes; Additional file 10).

Pan-GWAS analysis found a number of genes significantly associated with various phenotypic traits (Table 2). In terms of validating this analysis, one might naturally look for genes known to be associated with B. anthracis virulence [69] or B. cereus s. l.-induced food poisoning [22]; however, these genes are not found among the analysis results. Many of these genes were not annotated by Roary, and of the ones that were, some were not represented in the HAMSTR_FULL database, thus reducing the number of taxa for which there would have been usable data. The genes that were reported to be significantly associated with “disease: anthrax”, “disease: food poisoning”, and other traits thus represent hypotheses that remain to be validated. Only four traits had enough significant positively or negatively associated genes to allow for the identification of enriched subsets of genes involved in particular biological processes (“isolation source: cattle”, “isolation source: human”, “isolation source: non-primate mammal”, and “isolation source: soil”; Additional file 12). Of these, only the biological processes associated with “isolation source: soil” were sufficiently specific so as to be meaningfully interpretable. To increase the statistical power of the pan-GWAS analysis and thereby generate more comprehensive and specific lists of genes associated with various traits, one would need to include additional taxa with relevant metadata and gene content information.

All phylogenetic analyses in this study recapitulated the three-clade and seven-group classification systems, and taxa were consistently assigned to the same clade and group (Figs. 1, 2, 3 and 4), irrespective of the data source or analysis methodology used (Tables 3 and 4). This strongly suggests that the broad phylogenetic structure of B. cereus s. l. has been inferred correctly. I demonstrate that the three-clade and seven-group systems are compatible with each other, as no group has its member taxa assigned to multiple clades. Clades 1 and 2 are much more extensively sampled than Clade 3 due to historical interest in B. anthracis and B. thuringiensis (Additional file 8); a recent study has shown that there is likely to be a tremendous amount of as-yet incompletely characterized diversity in Clade 3 that can be assayed by sampling various natural environments [17]. Indeed, Clade 3 exhibited the greatest degree of species diversity; in particular, Group I contained representatives of seven different species, including two newly characterized species (B. bingmayongensis [6] and B. gaemokensis [7]; Additional file 8). Six of the 498 taxa did not place into one of the seven previously circumscribed groups, which suggests that classification systems will need to be updated and refined as additional isolates are sequenced. Perhaps most interesting among the unplaced taxa is B. manliponensis [8], which appears to be even more divergent from other B. cereus s. l. taxa than B. cytotoxicus [2] (Fig. 3 and Additional file 7). One possibility for updating the group-level classification system is to incorporate information from the Bayesian model-based clustering analysis, the results of which were shown to be compatible with the three-clade system and which recommended nine clusters instead of seven groups (Fig. 4).

Using the phylogeny of BCSL_498, I quantified the degree of monophyly for six current B. cereus s. l. species designations (Additional file 16 and Table 5). This analysis demonstrates quantitatively that with the exception of B. anthracis, species definitions within B. cereus s. l. are not currently based on phylogenetic relatedness, but rather on phenotypes such as virulence, physiology, and morphology [14, 18]. The primary focus of this study is the accurate reconstruction of phylogenetic relationships among taxa, and thus I make no specific recommendations for species re-designation based on these results. However, I do note a trend towards refined species designations that correlate with group affiliation; for example, several B. cereus strains in Group II have recently been re-designated B. wiedmanii [4]; similarly, Böhm et al. [18] suggested that all Group V taxa should be designated B. toyonensis [5]. In general, I recommend that taxonomic revisions are informed by well-supported phylogenetic hypotheses that have been generated without bias towards any particular species concept (e.g., dDDH boundaries [20]).

In a bioforensic setting, phylogenies that include well-supported strain-level relationships aid greatly in the identification of new isolates, and thus support both the attribution process (traceback of an isolate to its source) as well as analyses of pathogen evolution in an epidemic or outbreak scenario. However, the extremely high level of genomic conservation among closely related bacterial strains, especially in the core genome or in commonly typed conserved regions such as housekeeping genes, has limited the ability of previous analyses to make robust strain-level phylogenetic inferences. An important contribution of the current study is to show that bootstrap probabilities increase substantially when accessory genes are included in phylogenetic analyses along with core genes (Table 4). Thus, I have been able to resolve many strain-level, intra-group relationships of B. cereus s. l. with 100% bootstrap support for the first time (Additional file 14).

Conclusions

In this study, I used novel bioinformatic workflows to characterize the pan-genome and phylogeny of B. cereus sensu lato. Based on data from 114 complete genomes, I estimated that the B. cereus s. l. core and pan-genome contain ≈600 and ≈60,000 protein-coding genes, respectively. Pan-GWAS analysis revealed significant associations of particular genes with phenotypic traits shared by groups of taxa. All phylogenetic analyses recapitulated two previously used classification systems, and taxa were consistently assigned to the same major clade and group. By including accessory genes from the pan-genome in the phylogenetic analyses, I produced an exceptionally well-supported phylogeny of 114 complete B. cereus s. l. genomes. The best-performing methods were used to produce a phylogeny of all 498 publicly available B. cereus s. l. genomes, which was in turn used to compare three different classification systems and to test the monophyly status of various B. cereus s. l. species. The majority of the methodology used in this study is generic and could be leveraged to produce pan-genome estimates and similarly robust phylogenetic hypotheses for other bacterial groups. All scripts, software, and data products associated with this study are freely available.

Abbreviations

ACLAME: 

A CLAssification of mobile genetic elements

AFLP: 

Amplified fragment length polymorphism

BAPS: 

Bayesian analysis of population structure

BCSL: 

Bacillus cereus sensu lato

BLAST: 

Basic local alignment search tool

BP: 

Bootstrap probability

dDDH: 

digital DNA-DNA hybridization

GB: 

Gigabytes

gsi: 

Genealogical sorting index

GWAS: 

Genome-wide association study

HMM: 

Hidden Markov model

HaMStR: 

Hidden Markov model based search for orthologs using reciprocity

LCB: 

locally collinear block

MAFFT: 

Multiple alignment using fast fourier transform

MGE: 

Mobile genetic element

ML: 

Maximum likelihood

MLST: 

Multilocus sequence typing

MP: 

Maximum parsimony

NCBI: 

National center for biotechnology information

PANTHER: 

Protein ANalysis THrough evolutionary relationships

PATRIC: 

Pathosystems resource integration center

PAUP*: 

Phylogenetic analysis using parsimony *and other methods

PHYLIP: 

Phylogeny inference package

PRANK: 

Probabilistic alignment kit

RAM: 

Random access memory

RAxML: 

Randomized axelerated maximum likelihood

RF: 

Robinson-foulds

RefSeq: 

Reference sequence database

SNP: 

Single nucleotide polymorphism

Declarations

Acknowledgements

I thank Shashikala Ratnayake for assistance generating the Mash-distance-based phylogeny of Bacillus, Todd Treangen for helpful discussions about the project, and M.J. Rosovitz, Brian Janes, and Martina Eaton for providing feedback on drafts of the manuscript. I also thank two anonymous reviewers for their comments.

Funding

This work was funded under Contract No. HSHQDC-15-C-00064 awarded by the Department of Homeland Security (DHS) Science and Technology Directorate (S&T) for the operation and management of the National Biodefense Analysis and Countermeasures Center (NBACC), a Federally Funded Research and Development Center. The views and conclusions contained in this document are those of the author and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the DHS or S&T. In no event shall DHS, NBACC, S&T or Battelle National Biodefense Institute have any responsibility or liability for any use, misuse, inability to use, or reliance upon the information contained herein. DHS does not endorse any products or commercial services mentioned in this publication.

Availability of data and materials

All data analyzed during the current study were downloaded from public databases (ACLAME, NCBI, and PATRIC), and dates of download are provided in the text. A list of RefSeq assembly accessions for the taxa used in this study is provided in Additional file 1. HaMStR databases, concatenated phylogenetic data matrices, phylogenetic trees, Perl scripts and other data supporting the results of this study are available in the Dryad repository, Dryad DOI: doi:10.5061/dryad.dm82j.

Authors’ contributions

ALB designed and executed the study and wrote the manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The author declares that he has no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
National Biodefense Analysis and Countermeasures Center

References

  1. Nakamura LK. Bacillus pseudomycoides sp. nov. Int J Syst Evol Microbiol. 1998; 48(3):1031–5.Google Scholar
  2. Guinebretière MH, Auger S, Galleron N, Contzen M, De Sarrau B, De Buyser ML, Lamberet G, Fagerlund A, Granum PE, Lereclus D, De Vos P, Nguyen-The C, Sorokin A. Bacillus cytotoxicus sp. nov. is a novel thermotolerant species of the Bacillus cereus Group occasionally associated with food poisoning. Int J Syst Evol Microbiol. 2013; 63(1):31–40.View ArticlePubMedGoogle Scholar
  3. Lechner S, Mayr R, Francis KP, Prüss BM, Kaplan T, Wiessner-Gunkel E, Stewart GS, Scherer S. Bacillus weihenstephanensis sp. nov. is a new psychrotolerant species of the Bacillus cereus group. Int J Syst Evol Microbiol. 1998; 48(4):1373–82.Google Scholar
  4. Miller RA, Beno SM, Kent DJ, Carroll LM, Martin NH, Boor KJ, Kovac J. Bacillus wiedmannii sp. nov., a psychrotolerant and cytotoxic Bacillus cereus group species isolated from dairy foods and dairy environments. Int J Syst Evol Microbiol. 2016; 66(11):4744–53.View ArticlePubMedPubMed CentralGoogle Scholar
  5. Jiménez G, Urdiain M, Cifuentes A, López-López A, Blanch AR, Tamames J, Kampfer P, Kolsto A-B, Ramón D, Martínez JF, Codoner FM, Rosselló-Móra R. Description of Bacillus toyonensis sp. nov., a novel species of the Bacillus cereus group, and pairwise genome comparisons of the species of the group by means of ANI calculations. Syst Appl Microbiol. 2013; 36(6):383–91. doi:10.1016/j.syapm.2013.04.008.View ArticlePubMedGoogle Scholar
  6. Liu B, Liu GH, Hu GP, Cetin S, Lin NQ, Tang JY, Tang WQ, Lin YZ. Bacillus bingmayongensis sp. nov., isolated from the pit soil of Emperor Qin’s Terra-cotta warriors in China. Anton Leeuw. 2014; 105(3):501–10. doi:10.1007/s10482-013-0102-3.View ArticleGoogle Scholar
  7. Jung MY, Paek WK, Park IS, Han JR, Sin Y, Paek J, Rhee MS, Kim H, Song HS, Chang YH. Bacillus gaemokensis sp. nov., isolated from foreshore tidal flat sediment from the Yellow Sea. J Microbiol. 2010; 48(6):867–71. doi:10.1007/s12275-010-0148-0.View ArticlePubMedGoogle Scholar
  8. Jung MY, Kim JS, Paek WK, Lim J, Lee H, Kim PI, Ma JY, Kim W, Chang YH. Bacillus manliponensis sp. nov., a new member of the Bacillus cereus group isolated from foreshore tidal flat sediment. J Microbiol. 2011; 49(6):1027–32. doi:10.1007/s12275-011-1049-6.View ArticlePubMedGoogle Scholar
  9. Papazisi L, Rasko DA, Ratnayake S, Bock GR, Remortel BG, Appalla L, Liu J, Dracheva T, Braisted JC, Shallom S, Jarrahi B, Snesrud E, Ahn S, Sun Q, Rilstone J, Økstad OA, Kolstø A-B, Fleischmann RD, Peterson SN. Investigating the genome diversity of B. cereus and evolutionary aspects of B. anthracis emergence. Genomics. 2011; 98(1):26–39. doi:10.1016/j.ygeno.2011.03.008.View ArticlePubMedPubMed CentralGoogle Scholar
  10. Toby IT, Widmer J, Dyer DW. Divergence of protein-coding capacity and regulation in the Bacillus cereus sensu lato group. BMC Bioinformatics. 2014; 15(11):8. doi:10.1186/1471-2105-15-S11-S8.View ArticleGoogle Scholar
  11. Tettelin H, Masignani V, Cieslewicz MJ, Donati C, Medini D, Ward NL, Angiuoli SV, Crabtree J, Jones AL, Durkin AS, DeBoy RT, Davidsen TM, Mora M, Scarselli M, Margarit y Ros I, Peterson JD, Hauser CR, Sundaram JP, Nelson WC, Madupu R, Brinkac LM, Dodson RJ, Rosovitz MJ, Sullivan SA, Daugherty SC, Haft DH, Selengut J, Gwinn ML, Zhou L, Zafar N, Khouri H, Radune D, Dimitrov G, Watkins K, O’Connor KJB, Smith S, Utterback TR, White O, Rubens CE, Grandi G, Madoff LC, Kasper DL, Telford JL, Wessels MR, Rappuoli R, Fraser CM. Genome analysis of multiple pathogenic isolates of Streptococcus agalactiae: Implications for the microbial “pan-genome”. Proc Natl Acad Sci U S A. 2005; 102(39):13950–13955. doi:10.1073/pnas.0506758102. http://www.pnas.org/content/102/39/13950.full.pdf.View ArticlePubMedPubMed CentralGoogle Scholar
  12. Lapidus A, Goltsman E, Auger S, Galleron N, Ségurens B, Dossat C, Land ML, Broussolle V, Brillard J, Guinebretiere MH, Sanchis V, Nguen-the C, Lereclus D, Richardson P, Wincker P, Weissenbach J, Ehrlich SD, Sorokin A. Extending the Bacillus cereus group genomics to putative food-borne pathogens of different toxicity. Chem Biol Interact. 2008; 171(2):236–49. doi:10.1016/j.cbi.2007.03.003. Frontiers of Pharmacology and ToxicologyView ArticlePubMedGoogle Scholar
  13. Zwick ME, Joseph SJ, Didelot X, Chen PE, Bishop-Lilly KA, Stewart AC, Willner K, Nolan N, Lentz S, Thomason MK, Sozhamannan S, Mateczun AJ, Du L, Read TD. Genomic characterization of the Bacillus cereus sensu lato species: Backdrop to the evolution of Bacillus anthracis. Genome Res. 2012; 22(8):1512–24. doi:10.1101/gr.134437.111. http://genome.cshlp.org/content/22/8/1512.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  14. Guinebretière MH, Thompson FL, Sorokin A, Normand P, Dawyndt P, Ehling-Schulz M, Svensson B, Sanchis V, Nguyen-The C, Heyndrickx M, De Vos P. Ecological diversification in the Bacillus cereus group. Environ Microbiol. 2008; 10(4):851–65. doi:10.1111/j.1462-2920.2007.01495.x.View ArticlePubMedGoogle Scholar
  15. Tourasse NJ, Økstad OA, Kolstø A-B. HyperCAT: an extension of the SuperCAT database for global multi-scheme and multi-datatype phylogenetic analysis of the Bacillus cereus group population. Database. 2010; 2010:017. doi:10.1093/database/baq017.View ArticleGoogle Scholar
  16. Didelot X, Barker M, Falush D, Priest FG. Evolution of pathogenicity in the Bacillus cereus group. Syst Appl Microbiol. 2009; 32(2):81–90. doi:10.1016/j.syapm.2009.01.001.View ArticlePubMedGoogle Scholar
  17. Drewnowska JM, Swiecicka I. Eco-genetic structure of Bacillus cereus sensu lato populations from different environments in northeastern Poland. PLOS ONE. 2013; 8(12):1–11. doi:10.1371/journal.pone.0080175.View ArticleGoogle Scholar
  18. Böhm ME, Huptas C, Krey VM, Scherer S. Massive horizontal gene transfer, strictly vertical inheritance and ancient duplications differentially shape the evolution of Bacillus cereus enterotoxin operons hbl, cytk and nhe. BMC Evol Biol. 2015; 15(1):246. doi:10.1186/s12862-015-0529-4.View ArticlePubMedPubMed CentralGoogle Scholar
  19. Schmidt TR, Scott EJ, Dyer DW. Whole-genome phylogenies of the family Bacillaceae and expansion of the sigma factor gene family in the Bacillus cereus species-group. BMC Genomics. 2011; 12(1):430. doi:10.1186/1471-2164-12-430.View ArticlePubMedPubMed CentralGoogle Scholar
  20. Liu Y, Lai Q, Göker M, Meier-Kolthoff JP, Wang M, Sun Y, Wang L, Shao Z. Genomic insights into the taxonomic status of the Bacillus cereus group. Sci Rep. 2015; 5:14082.View ArticlePubMedPubMed CentralGoogle Scholar
  21. Okinaka RT, Keim P. The phylogeny of Bacillus cereus sensu lato. Microbiol Spectr. 2016; 4(1):TBS-0012-2012.Google Scholar
  22. Guinebretière MH, Velge P, Couvert O, Carlin F, Debuyser ML, Nguyen-The C. Ability of Bacillus cereus group strains to cause food poisoning varies according to phylogenetic affiliation (groups I to VII) rather than species affiliation. J Clin Microbiol. 2010; 48(9):3388–91. doi:10.1128/JCM.00921-10. http://jcm.asm.org/content/48/9/3388.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  23. O’Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, Rajput B, Robbertse B, Smith-White B, Ako-Adjei D, Astashyn A, Badretdin A, Bao Y, Blinkova O, Brover V, Chetvernin V, Choi J, Cox E, Ermolaeva O, Farrell CM, Goldfarb T, Gupta T, Haft D, Hatcher E, Hlavina W, Joardar VS, Kodali VK, Li W, Maglott D, Masterson P, McGarvey KM, Murphy MR, O’Neill K, Pujar S, Rangwala SH, Rausch D, Riddick LD, Schoch C, Shkeda A, Storz SS, Sun H, Thibaud-Nissen F, Tolstoy I, Tully RE, Vatsan AR, Wallin C, Webb D, Wu W, Landrum MJ, Kimchi A, Tatusova T, DiCuccio M, Kitts P, Murphy TD, Pruitt KD. Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 2015; 44(D1):733. doi:10.1093/nar/gkv1189.View ArticleGoogle Scholar
  24. Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, Phillippy AM. Mash: fast genome and metagenome distance estimation using MinHash. Genome Biol. 2016; 17(1):132. doi:10.1186/s13059-016-0997-x.View ArticlePubMedPubMed CentralGoogle Scholar
  25. Lefort V, Desper R, Gascuel O. FastME 2.0: A comprehensive, accurate, and fast distance-based phylogeny inference program. Mol Biol Evol. 2015; 32(10):2798–800. doi:10.1093/molbev/msv150. http://mbe.oxfordjournals.org/content/32/10/2798.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  26. Gascuel O. BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data. Mol Biol Evol. 1997; 14(7):685–95. http://mbe.oxfordjournals.org/content/14/7/685.full.pdf+html.View ArticlePubMedGoogle Scholar
  27. Ebersberger I, Strauss S, von Haeseler A. HaMStR: Profile hidden markov model based search for orthologs in ESTs. BMC Evol Biol. 2009; 9(1):157. doi:10.1186/1471-2148-9-157.View ArticlePubMedPubMed CentralGoogle Scholar
  28. Wattam AR, Abraham D, Dalay O, Disz TL, Driscoll T, Gabbard JL, Gillespie JJ, Gough R, Hix D, Kenyon R, Machi D, Mao C, Nordberg EK, Olson R, Overbeek R, Pusch GD, Shukla M, Schulman J, Stevens RL, Sullivan DE, Vonstein V, Warren A, Will R, Wilson MJC, Yoo HS, Zhang C, Zhang Y, Sobral BW. PATRIC, the bacterial bioinformatics database and analysis resource. Nucleic Acids Res. 2014; 42(D1):581–91. doi:10.1093/nar/gkt1099. http://nar.oxfordjournals.org/content/42/D1/D581.full.pdf+html.View ArticleGoogle Scholar
  29. Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014; 30(14):2068–9. doi:10.1093/bioinformatics/btu153. http://bioinformatics.oxfordjournals.org/content/30/14/2068.full.pdf+html.View ArticlePubMedGoogle Scholar
  30. Page AJ, Cummins CA, Hunt M, Wong VK, Reuter S, Holden MTG, Fookes M, Falush D, Keane JA, Parkhill J. Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics. 2015; 31(22):3691–3. doi:10.1093/bioinformatics/btv421. http://bioinformatics.oxfordjournals.org/content/31/22/3691.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  31. Löytynoja A, Goldman N. An algorithm for progressive multiple alignment of sequences with insertions. Proc Natl Acad Sci U S A. 2005; 102(30):10557–62. doi:10.1073/pnas.0409137102. http://www.pnas.org/content/102/30/10557.full.pdf.View ArticlePubMedPubMed CentralGoogle Scholar
  32. Price MN, Dehal PS, Arkin AP. FastTree 2 – approximately maximum-likelihood trees for large alignments. PLOS ONE. 2010; 5(3):1–10. doi:10.1371/journal.pone.0009490.View ArticleGoogle Scholar
  33. Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006; 23(2):254–67. doi:10.1093/molbev/msj030. http://mbe.oxfordjournals.org/content/23/2/254.full.pdf+html.View ArticlePubMedGoogle Scholar
  34. Huson DH, Steel M. Phylogenetic trees based on gene content. Bioinformatics. 2004; 20(13):2044–9. doi:10.1093/bioinformatics/bth198. http://bioinformatics.oxfordjournals.org/content/20/13/2044.full.pdf+html.View ArticlePubMedGoogle Scholar
  35. Bryant D, Moulton V. In: Guigó R, Gusfield D, (eds).NeighborNet: An Agglomerative Method for the Construction of Planar Phylogenetic Networks. Berlin: Springer; 2002, pp. 375–91. doi:10.1007/3-540-45784-4_28. http://dx.doi.org/10.1007/3-540-45784-4_28 Google Scholar
  36. Brynildsrud O, Bohlin J, Scheffer L, Eldholm V. Rapid scoring of genes in microbial pan-genome-wide association studies with Scoary. Genome Biol. 2016; 17(1):238. doi:10.1186/s13059-016-1108-8.View ArticlePubMedPubMed CentralGoogle Scholar
  37. Carbon S, Ireland A, Mungall CJ, Shu S, Marshall B, Lewis S, Hub A, Group WPW. AmiGO: online access to ontology and annotation data. Bioinformatics. 2008; 25(2):288. doi:10.1093/bioinformatics/btn615.View ArticlePubMedPubMed CentralGoogle Scholar
  38. Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, Thomas PD. PANTHER version 11: expanded annotation data from Gene Ontology and Reactome pathways, and data analysis tool enhancements. Nucleic Acids Res. 2016; 45(D1):183. doi:10.1093/nar/gkw1138.View ArticleGoogle Scholar
  39. Katoh K, Frith MC. Adding unaligned sequences into an existing alignment using MAFFT and LAST. Bioinformatics. 2012; 28(23):3144–6. doi:10.1093/bioinformatics/bts578. http://bioinformatics.oxfordjournals.org/content/28/23/3144.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  40. Eddy SR. Profile hidden Markov models. Bioinformatics. 1998; 14(9):755–63. doi:10.1093/bioinformatics/14.9.755. http://bioinformatics.oxfordjournals.org/content/14/9/755.full.pdf+html.View ArticlePubMedGoogle Scholar
  41. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool,. J Mol Biol. 1990; 215(3):403–10. doi:10.1006/jmbi.1990.9999.View ArticlePubMedGoogle Scholar
  42. Boto L. Horizontal gene transfer in evolution: facts and challenges. Proc R Soc Lond B Biol Sci. 2010; 277(1683):819–27. doi:10.1098/rspb.2009.1679. http://rspb.royalsocietypublishing.org/content/277/1683/819.full.pdf.View ArticleGoogle Scholar
  43. Coordinators NR. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2016; 44(Database issue):7–19. doi:10.1093/nar/gkv1290.Google Scholar
  44. Leplae R, Lima-Mendez G, Toussaint A. ACLAME: A CLAssification of Mobile genetic Elements, update 2010. Nucleic Acids Res. 2010; 38(suppl 1):57–61. doi:10.1093/nar/gkp938.View ArticleGoogle Scholar
  45. Birney E, Clamp M, Durbin R. GeneWise and Genomewise. Genome Res. 2004; 14(5):988–95. doi:10.1101/gr.1865504. http://genome.cshlp.org/content/14/5/988.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  46. Bazinet AL, Mitter KT, Davis DR, Van Nieukerken EJ, Cummings MP, Mitter C. Phylotranscriptomics resolves ancient divergences in the Lepidoptera. Syst Entomol. 2017; 42:305–16. doi:10.1111/syen.12217.View ArticleGoogle Scholar
  47. Bazinet AL, Cummings MP, Mitter KT, Mitter CW. Can RNA-Seq resolve the rapid radiation of advanced moths and butterflies (Hexapoda: Lepidoptera: Apoditrysia)? An exploratory study. PLOS ONE. 2013; 8(12). doi:10.1371/journal.pone.0082615.
  48. Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014; 30(9):1312–3. doi:10.1093/bioinformatics/btu033. http://bioinformatics.oxfordjournals.org/content/30/9/1312.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  49. Bazinet AL, Zwickl DJ, Cummings MP. A gateway for phylogenetic analysis powered by grid computing featuring GARLI 2.0. Syst Biol. 2014; 63(5):812–8. doi:10.1093/sysbio/syu031. http://sysbio.oxfordjournals.org/content/63/5/812.full.pdf+html.View ArticlePubMedPubMed CentralGoogle Scholar
  50. Pattengale ND, Alipour M, Bininda-Emonds ORP, Moret BME, Stamatakis A. In: Batzoglou S, (ed).How Many Bootstrap Replicates Are Necessary?. Berlin: Springer; 2009, pp. 184–200. doi:10.1007/978-3-642-02008-7_13. http://dx.doi.org/10.1007/978-3-642-02008-7_13 Google Scholar
  51. Sukumaran J, Holder MT. DendroPy: a Python library for phylogenetic computing. Bioinformatics. 2010; 26(12):1569–71. doi:10.1093/bioinformatics/btq228. http://bioinformatics.oxfordjournals.org/content/26/12/1569.full.pdf+html.View ArticlePubMedGoogle Scholar
  52. Bruen TC, Philippe H, Bryant D. A simple and robust statistical test for detecting the presence of recombination. Genetics. 2006; 172(4):2665–81. doi:10.1534/genetics.105.048975. http://www.genetics.org/content/172/4/2665.full.pdf.View ArticlePubMedPubMed CentralGoogle Scholar
  53. Treangen TJ, Ondov BD, Koren S, Phillippy AM. The Harvest suite for rapid core-genome alignment and visualization of thousands of intraspecific microbial genomes. Genome Biol. 2014; 15(11):524. doi:10.1186/s13059-014-0524-x.View ArticlePubMedPubMed CentralGoogle Scholar
  54. Didelot X, Wilson DJ. ClonalFrameML: Efficient inference of recombination in whole bacterial genomes. PLoS Comput Biol. 2015; 11(2):1–18. doi:10.1371/journal.pcbi.1004041.View ArticleGoogle Scholar
  55. Swofford DL. Phylogenetic analysis using parsimony (* and other methods). Version 4.Sunderland: Sinauer Associates; 2002.Google Scholar
  56. Robinson DF, Foulds LR. Comparison of phylogenetic trees. Math Biosci. 1981; 53(1):131–47. doi:10.1016/0025-5564(81)90043-2.View ArticleGoogle Scholar
  57. Rambaut A. http://tree.bio.ed.ac.uk/software/figtree/. Accessed Oct 2016.
  58. Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. 2016; 44(W1):242. doi:10.1093/nar/gkw290.View ArticleGoogle Scholar
  59. Cheng L, Connor TR, Sirén J, Aanensen DM, Corander J. Hierarchical and spatially explicit clustering of DNA sequences with BAPS software. Mol Biol Evol. 2013; 30(5):1224. doi:10.1093/molbev/mst028.View ArticlePubMedPubMed CentralGoogle Scholar
  60. Cummings MP, Neel MC, Shaw KL. A genealogical approach to quantifying lineage divergence. Evolution. 2008; 62(9):2411–22. doi:10.1111/j.1558-5646.2008.00442.x.View ArticlePubMedGoogle Scholar
  61. Bazinet AL, Cummings M. The Lattice Project: a Grid research and production environment combining multiple Grid computing models. Distributed & Grid Computing—Science Made Transparent for Everyone. Principles, Applications and Supporting Communities. Marburg: Tectum Publishing House; 2008, pp. 2–13.Google Scholar
  62. Fujie K, Hu HY, Tanaka H, Urano K, Saitou K, Katayama A. Analysis of respiratory quinones in soil for characterization of microbiota. Soil Sci Plant Nutr. 1998; 44(3):393–404. doi:10.1080/00380768.1998.10414461. http://dx.doi.org/10.1080/00380768.1998.10414461 View ArticleGoogle Scholar
  63. Kandeler E. Physiological and biochemical methods for studying soil biota and their function In: Stotzky G, editor. Soil Biochemistry: vol. 9. CRC Press, n.p.: 1996. p. 253–86. Chap. 7.Google Scholar
  64. Czaban J, Gajda A, Wróblewska B. The motility of bacteria from rhizosphere and different zones of winter wheat roots. Pol J Environ Stud. 2007; 16(2):301–8.Google Scholar
  65. Kaiser D. Bacterial swarming: A re-examination of cell-movement patterns. Curr Biol. 2007; 17(14):561–70. doi:10.1016/j.cub.2007.04.050.View ArticleGoogle Scholar
  66. Murphy SL, III RLT. In: Paul EA, (ed).Soil Microbiology, Ecology and Biochemistry: Academic Press, n.p.; 2006, pp. 53–83. Chap. 3.Google Scholar
  67. Bofkin L, Goldman N. Variation in evolutionary processes at different codon positions. Mol Biol Evol. 2006; 24(2):513. doi:10.1093/molbev/msl178.View ArticlePubMedGoogle Scholar
  68. World Health Organization and International Office of Epizootics and Food and Agriculture Organization of the United Nations. Anthrax in Humans and Animals. Nonserial Publication Series: World Health Organization, n.p.; 2008. https://books.google.com/books?id=EKYihvnaA7oC.
  69. Koehler TM. In: Koehler TM, (ed).Bacillus anthracis Genetics and Virulence Gene Regulation. Berlin: Springer; 2002, pp. 143–64. doi:10.1007/978-3-662-05767-4_7. http://dx.doi.org/10.1007/978-3-662-05767-4_7 Google Scholar

Copyright

© The Author(s) 2017

Advertisement