Complete rpoB gene sequences with homology to B. amyloliquefaciens DSM7^T were retrieved from the respective genomes of Bacillus strains available at NCBI (http://www.ncbi.nlm.nih.gov/sutils/genom_table.cgi?organism$=$microb). Sequence comparisons were obtained by NCBI BlastN (http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD$=$Web&PAGE_TYPE$=$BlastHome).

Alignment of DNA rpoB sequences was performed by the Clustal Omega program accessible at http://www.ebi.ac.uk/Tools/msa/clustalo/. A distance matrix was calculated from this alignment by DNA distance matrix calcuation (DNADIST program), and the matrix was then transformed into a tree by the NEIGHBOR program. In order to verify the accuracy of the tree multiple data sets were generated with the SEQBOOT program using 200 bootstrap replicates. A tree was built from each replicate with the DNADIST program, and then bootstrap values were computed with the CONSENSE program. The phylogenetic tree was visualized with TreeViewX (http://taxonomy.zoology.gla.ac.uk/rod/treeview.html). The programs used to construct the phylogenetic tree were obtained from the PHYLIP package, v.3.65 (Felsenstein, 1989), which is accessible at http://evolution.genetics.washington.edu/phylip.html.

Comparative genome analysis was performed using the EDGAR 1.3 software framework. For orthology estimation EDGAR uses a generic orthology threshold calculated from the similarity statistics of the compared genomes (Blom et al., 2016; http://edgar.computational.bio.uni-giessen.de). A private project was constructed comprising 66 genomes closely related to B. amyloliquefaciens DSM7^T and selected other representatives of the B. subtilis species complex. To construct a phylogenetic tree for this project, around 2000 core genes were computed by pairwise iterative comparison of a set of genomes (Blom et al., 2016). In a following step multiple alignments of the core genes were generated using MUSCLE, non-matching parts of the alignment were masked by GBLOCKS and subsequently removed. The remaining parts of all alignments were concatenated to one large alignment. The PHYLIP package was used to generate a phylogenetic tree of this alignment, represented in newick format.

The EDGAR software framework was also used to calculate average nucleotide identity (ANI) and average amino acid identity (AAI), matrices for a selected set of genomes. The blast hits between the orthologous genes of the core of the selected genome were analyzed for their mean/median percent identity values. The recommended species cut-off was 95% for the ANI and AAI indices (Richter and Rosselló-Móra, 2009). In addition, JSpeciesWS (http://jspecies.ribohost.com/jspeciesws/) was used to determine ANIb (average nucleotide identity based on BLAST+) and ANIm (average nucleotide identity based on MUMmer) values by pairwise genome comparisons. Correlation indexes of their Tetra-nucleotide signatures (TETRA) were determined by using the JSpeciesWS software (Richter et al., 2016).

The genome-to-genome-distance calculator (GGDC) version 2.1 provided by DSMZ (http://ggdc.dsmz.de/) was used for genome-based species delineation (Meier-Kolthoff et al., 2013) and genome-based subspecies delineation (Meier-Kolthoff et al., 2014). Distances were calculated by (i) comparing two genomes using the chosen program to obtain HSPs/MUMs and (ii) inferring distances from the set of HSPs/MUMs using three distinct formulas. Next, the distances were transformed to values analogous to DDH. The DDH estimates were based on an empirical reference dataset comprising real DDH values and genome sequences. The DDH estimate resulted from a generalized linear model (GLM) which also provided the estimate's confidence interval (after the ± sign). Three formulas are available for the calculation: Formula: 1 (HSP length/total length), formula: 2 (identities/HSP length) and formula 3 (identities/total length). Formula 2, which is especially appropriate to analyze draft genomes, was used.

The core genomes of 20 type strains of the B. subtilis species complex were used for phylogenomic analysis applying the EDGAR software package (Figure 1). Four main monophyletic groups were corroborated by 100% bootstrap values. Clade I (“subtilis”) is early diverged into two branches comprising B. atrophaeus, and B. subtilis and its close relatives; clade II (“amyloliquefaciens”) comprises B. amyloliquefaciens, B. siamensis, and a conspecific group containing the type strains of B. amyloliquefaciens subsp. plantarum, B. velezensis, and B. methylotrophicus; clade III (“licheniformis”) consists of B. licheniformis and B. sonorensis; and clade IV (“pumilus”) comprises B. pumilus, B. safensis, B. xiamenensis, and a conspecific group involving the type strains of B. altitudinis, B. stratosphericus, and B. aerophilus. The members of clade II appeared closely related. This is indicated by the high number of orthologous CDSs (2794) shared by the five type strains of clade II. A similar cladogram has been published recently (Dunlap C. et al., 2016) suggesting that the B. subtilis species complex can be divided into four groups above species level, which need further characterization. We have directed our further analysis to clade II (named from now on “operational group B. amyloliquefaciens”), which clearly shows the highest degree of compactness.

It is obvious, that 16S rRNA sequences are not sufficient to discriminate representatives of the B. subtilis species complex. For example, comparison of the complete 16S rRNA sequences of B. amyloliquefaciens DSM7^T and B. subtilis 168^T revealed 99.48% identity (Table 1), which is well above of the recommended threshold of >98.65% for species delineation (Kim et al., 2014). In order to elucidate more precisely the phylogenetic and taxonomic relationship of the members of the B. subtilis species complex belonging to the “operational group B. amyloliquefaciens,” we used two methods. (i) Tetra correlation search (TCS, Richter et al., 2016) was performed with the complete genome of DSM7^T and (ii) the complete RNA polymerase beta-subunit (rpoB) gene of DSM7^T was used for BLASTN comparison with the corresponding sequences extracted from complete genomes or genome assemblies. Fifty-Two genomes, which were in range with the intraspecific Tetra-nucleotide signature correlation index (>0.99) were detected in the JSpecies data bank. The TCS value determined for B. subtilis was only 0.954, suggesting that using this alignment-free parameter allows discriminating of B. subtilis and B. amyloliquefaciens (Table 1). Complete rpoB gene sequencing has been proposed as phylogenetic marker (Klenk et al., 1994) and as a supplement to DDH (Adékambi et al., 2008). The power and potential of complete rpoB gene sequence in taxonomic, phylogenetic and evolutionary studies has been previously reported (Sharma and Patil, 2011). Our BLASTN search revealed that at least 66 genomes present in the NCBI data bank contain rpoB gene sequences with more than 98% identity to the rpoB gene from DSM7^T, the type strain of B. amyloliquefaciens (Priest et al., 1987). For comparison, the rpoB gene from B. subtilis subsp. subtilis 168^T displayed only 90.3% identity to B. amyloliquefaciens. The rpoB gene identities among strains assigned as being B. amyloliquefaciens, B. siamensis, B. amyloliquefaciens subsp. plantarum, B. methylotrophicus, B. velezensis, and B. vanillea are listed in Table 1. The list of strains containing rpoB genes with high similarity to B. amyloliquefaciens DSM7^T includes also strains obviously not correctly assigned, such as B. subtilis, Bacillus sp., or Paenibacillus polymyxa. It is interesting to note that majority of the strains representing the conspecific B. velezensis/B.methylotrophicus/B.amyloliquefaciens subsp. plantarum group were isolated from plant sources, whilst B. amyloliquefaciens sensu stricto seems to be soil-borne. The main source of the salt tolerant B. siamensis/B.vanillea group was fermented plant food (Table 1).

The phylogenetic tree based on complete rpoB gene sequence suggests existence of three tightly connected monophyletic groups: (i) B. amyloliquefaciens containing six strains including type strain DSM7^T; (ii) B. siamensis cluster consists of three strains: the type strain KCTC 13613^T, strain XY18, originally assigned as type strain for B. vanillea (Chen et al., 2014) but recently reclassified as being B. siamensis (Dunlap, 2015), and a strain assigned as being B. amyloliquefaciens JJC33M; (3) the conspecific complex comprising B. velezensis, B. methylotrophicus, and B. amyloliquefaciens subsp. Plantarum contained 57 strains. The tree is robust displaying high bootstrap values for all three groupings, although the three clusters are closely related and separated by only 0.01–0.02 substitutions per nucleotide position. By contrast, taxonomic distance to B. subtilis is around tenfold larger (Figure 2).

In order to confirm the phylogenetic analysis based on rpoB sequences we calculated the core genomes using the EDGAR 1.3 program package. A total of 1998 CDSs were shared by the 66 core genome sequences extracted in that analysis. It ruled out that the phylogenomic tree based on complete core genome sequences (Figure 3) did reflect the phylogenomic distances similar as the phylogenetic tree based on rpoB nucleotide sequences (Figure 2). The same robust monophyletic groups as in Figure 2 were obtained. The B. siamensis cluster consisting of three representatives shared a core genome of 3097 CDSs; the B. amyloliquefaciens cluster consisting of six representatives shared a core genome of 3139 CDSs; and the conspecific group containing 57 plant-growth promoting Bacilli including FZB42^T shared a relatively small core genome consisting of only 2295 CDSs, which is mainly due to the high number of genomes included in this analysis. Subgroups of this cluster shared core genomes ranging from 2659 to 3137 CDSs (Figure 3). Again, the NJ tree suggested that B. amyloliquefaciens subsp. Plantarum, B. methylotrophicus, and B. velezensis formed a monophyletic group corroborating recent findings (Dunlap C. A. et al., 2015; Wu et al., 2015; Dunlap C. et al., 2016).

At next we tried to elucidate the taxonomic status of these closely related genomes. Different phylogenetic and phylogenomic methods were used to analyze relationship of all 65 genomes with that of B. amyloliquefaciens DSM7^T. As shown above, rpoB sequence similarity, exceeding threshold of species delineation, and the intraspecific Tetra-nucleotide signature correlation index (>0.99) suggested that all strains analyzed belong to the species B. amyloliquefaciens. TETRA analysis (Jspecies) demonstrated that the six type strains of clade II were closely related and yielded pairwise Tetra results (tetranucleotide signature correlation index) in species range (≥0.989, Figure 4 lower part). Deviations of the mean G+C content calculated for the whole genomes were less than one percent which does not contradict species definition (Table 2). Grouping of all strains into a single species, B. amyloliquefaciens, was further supported by the AAI values (Table 1). The mean AAI values of the 66 core genomes selected by their rpoB similarity to DSM7^T were ≥96.5%, exceeding the proposed cut-off of 96% for species delineation. However, parameters, considered recently as being most important for genome-based species delineation, such as ANI and dDDH (Federhen et al., 2016), did not support this conclusion (Table 1).

ANI analysis performed with the EDGAR program package discriminated clearly two clusters corresponding to clades I and II of the B. subtilis species complex (Figure 4). Clade II representing the B. amyloliquefaciens group was divided into three groups consisting of B. amyloliquefaciens DSM7^T (i), B. siamensis and B. vanillea (ii), and the conspecific complex formed by the type strains B. methylotrophicus KACC13105^T, B. velezensis KCTC 13102^T and B. amyloliquefaciens subsp. plantarum FZB42^T (iii). The latter group displayed ANI values of >98% exceeding the cut-off for species delineation when compared with each other suggesting that the members of the conspecific complex belong to a single species. B. methylotrophicus KACC 13105^T, B. velezensis KCTC 13102^T, and B. amyloliquefaciens subsp. plantarum FZB42^T displayed similar median ANI values ranging between 94.3 and 94.8% when compared with B. amyloliquefaciens DSM7^T (1) and B. siamensis KCTC-13613^T (2), respectively. Given a calculated deviation of ±2.2–2.3% the ANI matrix values suggests a high degree of relatedness to B. amyloliquefaciens, B. siamensis and the conspecific group formed by B. amyloliquefaciens subsp. plantarum, B. methylotrophicus, and B. velezensis, but did not sufficiently support species delineation (Figure 4, upper part). According to more recent findings the recommended cut-off point for species delineation corresponds to ~96% ANI (Colston et al., 2014). Similar results were obtained when ANIb and ANIm values were determined by using the JSpecies program package for all the 66 genomes included in this analysis. Threshold values sufficient for species delineation were only obtained, when representatives of B. amyloliquefaciens (6 genomes), B. siamensis (3 genomes), and of the conspecific group (57 genomes) were compared with their respective type strains. However, comparison of the 57 strains of the conspecific group (e.g., FZB42, B. velezensis) with either B. amyloliquefaciens DSM7^T or B. siamensis KCTC 13613^T yielded ANI values slightly below the cut-off for species delineation. The same was true when the three members of the B. siamensis group were compared with either FZB42^T or DSM7^T (Table 2) suggesting that according to ANI analysis the members of clade II represent three discrete, although closely related, species.

In order to finally decide, whether all strains of clade II belong to one species or not, electronic DNA-DNA hybridization (dDDH) was applied in a quantitative analysis involving all 66 genomes. As shown previously, dDDH is useful to mimic the wet-lab DDH and can be used for genome-based species delineation and genome-based subspecies delineation (Meier-Kolthoff et al., 2013, 2014). For calculating dDDH three different formulas can be applied (see Materials and Methods), but only results obtained with the recommended formula 2 were used in our analysis (Table 2). When comparing members of the “siamensis group 2” and the “conspecific B. velezensis group” with B. amyloliquefaciens DSM7^T, dDDH values of <70%, the defined threshold for species delineation, were obtained. All in all, dDDH supports our previous finding about a close relationship within clade II, but did not support their classification into one single species. The results are summarized in Table 2 and Supplementary Table 1.

Compared to other members of the B. subtilis species complex, the plant-associated B. amyloliquefaciens possess an enormous potential to synthesize bioactive secondary metabolites. Besides five gene clusters, known from B. subtilis to mediate nonribosomal synthesis of secondary metabolites, four giant gene clusters absent in B. subtilis 168 were identified in FZB42 (Chen et al., 2007). The nine gene clusters that direct the synthesis of bioactive peptides and polyketides by modularly organized mega-enzymes define both nonribosomal peptide synthetases (NRPSs) and polyketide synthases (PKS). Three (bmyD, dfn, and mln) are not present in B. subtilis 168, but occur in all members of the “operational group B. amyloliquefaciens.” Except for the gene cluster encoding bacilysin synthesis, the functional activities of the remaining gene clusters depend on Sfp, an enzyme that transfers 4′-phosphopantetheine from coenzyme A to the carrier proteins of nascent peptide or polyketide chains. A direct comparison revealed that the nine gene cluster responsible for nonribosomal synthesis of bioactive secondary metabolites including macrolactin are only present in FZB42 and in the other members of the conspecific B. velezensis group, whilst the gene cluster involved in macrolactin synthesis was not detected in B. siamensis and B. amyloliquefaciens (Table 3). Noteworthy, the gene cluster responsible for synthesis of the polyketide difficidin was present in B. siamensis, but not in any other member of the B. subtilis species complex suggesting a stepwise loss of the ability to synthesize secondary metabolites in the order B. velezensis (including FZB42) B. siamensis B. amyloliquefaciens.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.