Bacillus tequilensis EA-CB0015 was deposited in ATCC as PTA-123533. The strain was revived from glycerol stocks onto half strength tryptic soy agar [20 g/L of TSA (Oxoid) and 9 g/l of Bacto agar (BD)] and incubated at 37°C for 24 h. A single colony was inoculated into a 100 mL Erlenmeyer containing 20 mL of tryptic soy broth (Oxoid) and grown for 8 h in a shaking incubator (250 rpm, 37°C). Genomic DNA was extracted using the UltraClean Microbial DNA Isolation kit (MoBio) for Illumina sequencing and DNeasy Blood & Tissue kit (Qiagen) for Single Molecule Real-Time (SMRT) sequencing. Samples were sequenced on an Illumina HiSeq 2000 (100 bp paired end reads) and PacBio RS II (Macrogen, Inc). Contigs less than 100 bp were discarded. Unicycler v0.4.8.0 was used for hybrid assembly using the normal bridging mode (Wick et al., 2017).

The assembled genome was annotated using the NCBI Prokaryotic Genome Annotation Pipeline (PGAP) v6.0 (Haft et al., 2018) and deposited under accession no. NZ_CP048852. SkewDB was used to predict replication origins and termini (Hubert, 2022). Genome figures were created using the CGView Server (Grant and Stothard, 2008). Predicted proteins were classified into clusters of orthologous groups (COGs) using WebMGA (Wu et al., 2011), and analyzed in the Kyoto Encyclopedia of Genes and Genomes database (KEGG; entry T07089) to assign roles in metabolic pathways (Kanehisa et al., 2016).

All strains analyzed in this study are listed in Supplementary Table S1. Reference strains included B. subtilis species complex (B. pumilus, B. subtilis, B. tequilensis, B. mojavensis, B. atrophaeus, B. velezensis, and B. amyloliquefaciens, B. licheniformis), the B. cereus species complex (B. cereus, B. thuringiensis, B. anthracis, B. mycoides), and B. coagulans. Clostridium kluyveri DSM 555^T served as the outgroup.

Genomes were downloaded from public sequence repositories and annotated using Prokka v1.14.6 (Seemann, 2014). Genes encoding groEL (chaperonin, large subunit), gyrA (DNA gyrase, subunit A), rpoB (RNA polymerase, beta subunit), polC (DNA polymerase III, alpha subunit), and purH (phosphoribosylaminoinidazole carboxamide formyltransferase) were retrieved from each strain, aligned using MAFFT v7.4 (Katoh, 2002), manually trimmed, and concatenated head to tail using AMAS v0.98 (Borowiec, 2016). The matrix (12,483 positions; 1,647 from groEL, 2,517 from gyrA, 3,030 from polC, 1,542 from purH, and 3,747 from rpoB) was analyzed using ModelTest-NG v0.1.6 (Darriba et al., 2020) to determine the best substitution model for each partition. Phylogenetic trees were constructed using maximum-likelihood (ML) and Bayesian inference (BI) methods. The ML phylogeny using the GTR + I + G4 model was best for groEL, gyrA, purH and rpoB, and GTR + G4 for polC, with bootstrapping (1,000) using IQ-TREE v1.6.9 (Nguyen et al., 2015). The BI phylogeny was calculated using Markov Chain Monte Carlo (MCMC) analysis and the GTR + I + G4 model using Mr. Bayes v.3.2 (Ronquist et al., 2012). Two independent runs were performed for 10 million generations (1 sampling every 1,000 generations). Effective sample size, convergency, and stationarity values were evaluated using Tracer v1.7.1 (Rambaut et al., 2018). The final phylogenetic tree was generated using DendroPy v4.4.0 (Sukumaran and Holder, 2010) and visualized using iToL v5.0 (Letunic and Bork, 2019).

The ability of strain EA-CB0015 to utilize different carbon sources was performed using a BioLog GEN III Microplate (Biolog Inc.). EA-CB0015 was grown on Universal Growth Agar (BUG; Biolog Inc.) and incubated at 30°C for 24 h. Isolated colonies were resuspended in the “Inoculation Fluid-A” (IF-A; Biolog Inc.) following manufacturer’s instructions and this suspension was inoculated into the Gen III microplate. After 24 h incubation at 30°C, the microplate was analyzed qualitatively for color development using a MicroStation™ 2 Reader (Biolog Inc.), measuring OD at 595 and 750 nm.

Colonies grown on TSA plates were stained for endospores using Schaeffer-Fulton method (Gerhardt et al., 1994). Micrographs were recorded using an Axiostar Plus microscope (Carl Zeiss) at 100X equipped with an AxioCam ICc3 (Carl Zeiss) and ZEN 2.3 lite software (Carl Zeiss). Biofilm formation was examined by inoculating a diluted culture of EA-CB0015 into 12-well plates with LBGM media. These were incubated for 48 h at 30°C (statically) and photographed (7.1 MP digital camera, Kodak EasyShare P712). Motility was assayed by inoculating 3 μl EA-CB0015 grown on tryptic soy broth TSB (Oxoid; OD 1.0, 600 nm) onto the center of 10% LB plates solidified with either 0.7% (swarming) or 0.3% agar (swimming), and incubated at 30°C for 24 h (Kearns and Losick, 2003; Ghelardi et al., 2012). Flagella were visualized using Ryu stain (Heimbrook et al., 1989).

Genomes were analyzed using antiSMASH v5.1.2 (Blin et al., 2019). BGCs with similarity scores > 70% were reported for all 58 Bacillus strains. Regions encoding with similarity scores less than <70% were further analyzed for strain EA-CB0015 using BLAST (Johnson et al., 2008) and HMMER v3.3.1 (Eddy, 2011). BGCs for natural products known to be produced by Bacilli but not detected by antiSMASH were manually identified by BLAST.

Putative phage and phage-like regions were identified using PHASTER (Zhou et al., 2011; Arndt et al., 2016). Data ordering, sub-setting and reshaping was performed using Tidyverse v1.3.0 (Wickham et al., 2019) and an in-house script uploaded on GitHub¹ (García-Botero et al., 2021). GEIs were identified using IslandViewer 4 (Bertelli et al., 2017) and ISs identified using ISfinder (e-value cutoff of 1e-07; Siguier et al., 2006). Integrative or conjugative elements (ICEBs1) were predicted with ICEBerg 2.0 (Liu M. et al., 2019). Antimicrobial resistance genes were identified using ResFinder-4.0 (70% identity threshold; minimum length of 70%; Bortolaia et al., 2020), and putative virulence factors from the Virulence Factors of Pathogenic Database (VFDB; Liu B. et al., 2019). RM systems were retrieved from REBASE (Roberts et al., 2015). Cluster BLAST analysis of identified RM systems were performed using Cblaster v1.3.0 (Gilchrist et al., 2021). Potential TA systems identified using TAfinder (Xie et al., 2018) and manually BLAST searches. Finally, CRISPRCasFinder (Couvin et al., 2018) was used to identify CRISPRs and Cas genes.

Synteny was analyzed using Easyfig v2.2.2 (Sullivan et al., 2011). ANI values (Richter and Rosselló-Móra, 2009) were calculated using GTDB-Tk v0.3.2 (Chaumeil et al., 2019). Comparisons between B. subtilis 168, B. tequilensis ATCC BAA 819^T and B. tequilensis EA-CB0015 were performed using proteome and protein family functions in PATRIC (Davis et al., 2020). Presence or absence of predicted genes and proteins in EA-CB0015 were compared to B. subtilis 168 (NBCI Reference Sequence: NC_000964.3) using SubtiWiki (Zhu and Stülke, 2018). E-values < 1e-5 and sequence identities > 70% were used as cutoffs.

Hybrid assembly of Illumina and SMRT sequencing reads produced a single circularized chromosome of 4,012,371 bp with an average GC content of 43.7%. In total, 4,112 CDS (3,951 genes, 161 pseudogenes), 10 copies of rRNAs (5S, 16S, and 23S), 86 tRNAs, and 5 ncRNAs were present. The forward strand encoded 2,134 genes and the reverse strand 1,978 genes (Figure 1). GC-skew suggests single symmetrical, bi-directional replication of the genome (Supplementary Figure S1). COG analysis assigned these into 5,163 protein families, most of which belonged to general or unknown function (categories R and S). Besides these, the largest groups of proteins are dedicated for amino acid transport and metabolism (E), energy production and conversion (C), translation, ribosomal structure, and biogenesis (J), DNA replication, recombination, and repair (L), and transcription (K; Supplementary Figure S2).

EA-CB0015 was originally identified as B. subtilis (Ceballos et al., 2012) and later referred to as B. tequilensis (Cuellar-Gaviria et al., 2021). This ambiguity derived from the use of 16S rRNA gene phylogenies, which do not provide sufficient resolution at the species level for strains within the B. subtilis complex (Rooney et al., 2009). To definitively determine the species identity of the strain, we performed multilocus sequence analysis (MLSA) using established marker genes shown to effectively determine speciation in Bacillus. EA-CB0015 was more closely related to B. tequilensis ATCC BAA 819^T than to other B. subtilis subspecies, including subtilis, natto, stercoris, spizizenii and inaquosorum (Figure 2A). This is supported by high ANI values between the genomes of strain EA-CB0015 and B. tequilensis ATCC BAA 819^T (98.6%), vs. B. subtilis subsp. inaquosorum KCTC 13429^T (92.3%) B. subtilis subsp. subtilis 168 (91.4%; Supplementary Table S2). Thus, EA-CB0015 belongs to the phylogenetically homogeneous B. subtilis species complex and is classified as B. tequilensis.

The phylogenetic analysis suggests B. tequilensis and B. subtilis share a close evolutionary history. To understand the nature of this relationship, we identified genetic similarities and differences between EA-CB0015, ATCC BAA 819^T, and 168. The GC content of all three strains were similar (43.5%–44.0%) and they contained the identical number of rRNA and tRNA encoding genes. The B. tequilensis genomes were comparable in size (4 Mbp) and the number of protein coding sequences, but both were smaller than B. subtilis 168 (4.2 Mbp; Supplementary Table S3). Except for an extra 1.4 kbp plasmid in ATCC BAA 819^T, both B. tequilensis genomes were largely syntenic (Figure 3A). Although less conservation was observed between 168 and EA-CB0015, large segments of their genomes remained homologous with each other. A 132 kbp inversion corresponding to bacteriophage SPβ was near the replication terminus (terC) of both strains (Figure 3A). In strain 168, the prophage was encoded on the complementary strand 134 kb upstream terC, while in strain EA-CB0015 it was on the positive strand 186 kb downstream of terC. Interestingly, this region was highly variable between 168 and EA-CB0015, and completely absent from the genome of B. tequilensis ATCC BAA 819^T. Compared to 168, the SPβ prophage region of EA-CB0015 lacked the sublancin BGC and its immunity gene (Denham et al., 2019), a gene encoding for DNA cytosine-5-methyltransferase (mtbP), and numerous hypothetical proteins with unknown function. Most of the encoded proteins and their families were shared between all three strains (Figure 3B; Supplementary Table S4).

Multicellular behavior is used by B. subtilis to afford survival in its natural environment. These mechanisms include motility, biofilm formation, cannibalism, quorum sensing, competence, and sporulation (González-Pastor, 2017). As B. subtilis 168 is well studied for its ability to assemble multicellular communities and closely related to EA-CB0015 (O’Toole et al., 2000; González-Pastor, 2017), we examined if genes for these phenotypes were also present in our strain.

Genes for flagellum-mediated motility including chemotaxis and swarming were present in the EA-CB0015 genome. All thirty-two genes required to form the flagellar basal body and hook were encoded within a 27-kb fla/che operon, as well as all the genes needed for filament formation (flagellin monomer protein Hag, FliD, FlgK, FlgL), assembly, and rotation (Supplementary Table S5A; Mukherjee and Kearns, 2014). Other genes required for the basal body formation (FlhO and FlhP), torque in flagellar rotation (MotA, MotB), cell separation (LytC), bistable regulation (SigD, SwrA, SwrB, DegS/DegU, YmdB, and SlrA/SinR/SlrR), chemotaxis and associated chemoreceptors (2 soluble and 7 membrane-bound) were also present (Supplementary Tables S5A,B). Microscopy showed motile, polarly flagellated cells (Supplementary Figure S3). These results suggest EA-CB0015 can sense and swim toward a wide range of attractants (pH, amino acids, oxygen, etc.) that are present in its local environment.

To form a biofilm, cells must transition from motile to sessile states, aggregate, and embed themselves within a self-produced extracellular matrix (EM; Arnaouteli et al., 2021). EA-CB0015 formed biofilms in static culture (Supplementary Figure S3), and as expected, the genome of EA-CB0015 encoded genes essential for this process including exopolysaccharide biosynthesis (epsA-O), protein fiber TasA (tapAsipWtasA operon), the hydrophobin-like protein BslA, and genes for poly-γ-glutamic acid (γ-PGA) biosynthesis (pgsBCAE; Supplementary Table S5C). Production of EM is also linked to multiple regulatory proteins that were all present within the genome (Supplementary Table S5C). The Spo0A pathway (KinC, KinD, Spo0A, SinI/SinR/SlrR, AbrB) controls the expression of more than 100 genes including biofilm matrix gene expression and sporulation, while the YwcC-SlrA stress response pathway facilitates adaptation to changing environmental conditions. The DegS–DegU two-component system regulates competence, motility and secretion of degradative enzymes. Lastly, genes that mediate expression of slrR, an essential regulator of biofilm formation were also present (abh, ymdB, remA, remB; Vlamakis et al., 2014).

Interestingly, EA-CB0015 does not harbor the yitPOM operon. The paralogous sdpABC operon, is also absent from the genome of the strain. The yitPOM operon expresses the biofilm associated toxin (YitM) and the extracellular protease (NprB), both of which suppress competitors in B. subtilis biofilms (Kobayashi and Ikemoto, 2019). In addition, surfactin, encoded in the srfA operon (Supplementary Table S5B), may induce potassium leakage that stimulates the sensor kinase KinC. This may also activate expression of biofilm formation genes (López et al., 2009).

Cannibalistic behavior delays the entry into sporulation within a subpopulation of cells. This is controlled by the production of sporulation delaying protein-Sdp and sporulation killing factor-Skf, which lyse and kill sensitive siblings (González-Pastor et al., 2003). Interestingly, EA-CB0015 contained genes for Skf (skfA-H), but not Sdp (sdpABC and sdpRI operons were absent; Supplementary Table S5D). These results suggest EA-CB0015 cells may exhibit an accelerated sporulating phenotype. As Spo0A-inactive cells are not lysed, the pool of nutrients released into the environment from cell death is reduced. Thus, neighboring cells experience starvation to initiate spore formation (González-Pastor et al., 2003).

Genes encoding mechanisms of quorum sensing were present in the EA-CB0015 genome (Supplementary Table S5E). These may allow strain EA-CB0015 to coordinate physiological processes such as the synthesis of exoproteases and other extracellular enzymes in response to cell density. The lipopeptide surfactin is also positively regulated by the phosphorylated form of ComA, which is part of the quorum sensing (QS) system (ComQXPA) in B. subtilis (Nakano et al., 1991; Kalamara et al., 2018). Interestingly, the putative ComQ, ComX, and ComP proteins shared low sequence identity with their corresponding homologs in 168 (Supplementary Table S5E). This suggests the Com system may handle separate social communication groups or pherotypes (Stefanic and Mandic-Mulec, 2009; Oslizlo et al., 2014). The genome also encoded six putative receptor-signal pairs of the Rap-Phr system (Rap-Phr A, C, E, F, H; Supplementary Table S5E). Similar to 168, accumulation of Phr peptides in EA-CB0015 may suppress effects of Rap proteins to allow expression of genes for swarming motility, biofilm formation, exoprotease secretion and genetic competence (Kalamara et al., 2018). Related to quorum sensing, EA-CB0015 also contained genes for the acquisition and incorporation of extracellular DNA into the host cell (genetic competence). These include the master competence regulator comK, genes encoding proteins essential for DNA binding and import (comC, operons comE, comF, comG, bdbD, bdbC), and cytosolic proteins that modulate recombination and transformation efficiency (RecA, SsbB, DprA, CoiA, NucA). Genes encoding proteins for transcriptional (Rok, CodY, Kre) and post-translational regulation (MecA, ComS) of ComK were also present (Supplementary Table S5F).

Endospores were visible when cells were subjected to Ryu stain (Supplementary Figure S3). Indeed, genes needed for sporulation were present in the genome (Supplementary Table S5G). These included genes encoding for morphogenetic proteins (SpoIVA, SpoVM, SpoVID, SafA, CotE, CotX/CotY/CotZ), their interacting partners during spore coat assembly, and spore crust proteins. Interestingly, some shared < 60% sequence identity with the homologs in 168, while others were completely absent (Supplementary Table S5G). The biosynthetic genes for legionaminic acid, which is used in crust formation, were notably missing (Supplementary Table S5G). Their absence suggests that the surface of EA-CB0015 spores may be reduced in hydrophilicity and charge (Dubois et al., 2020).

Most of the genes encoding signal transduction proteins (histidine sensor kinases KinA-E, master regulator Spo0A, phosphotransferases Spo0F and Spo0B, etc) related to sporulation were identified (Supplementary Table S5G). Homologs of lrpAB were missing, but other studies have indicated negligible effects on glyA transcription or sporulation through KinB (Dartois et al., 1997). The absence of sivC may result in a greater sporulation efficiency, as it functions as an inhibitor of the KinB and KinC pathway (Garti-levi et al., 2013). As expected, genes encoding proteins involved with major events in spore gemination were also present (Setlow et al., 2017). These include genes needed for germinant sensing (germinant receptors), release of dipicolinic acid (DPA) related to heat resistance (spoVA proteins, GerD), and hydrolysis of cortex peptidoglycan (CwlJ, SleB, SleL; Supplementary Table S5G). Altogether, these data suggest EA-CB0015 endospores can fully germinate once in favorable environmental conditions.

Nutrient availability on leaves, especially for organic compounds, is spatially heterogeneous and limited (Leveau and Lindow, 2001; Lindow and Brandl, 2003). A major carbon source that leach from the interior of the plant are common sugars (Wildman and Parkinson, 1981; Fiala et al., 1990; Lindow and Brandl, 2003). EA-CB0015 encodes all the genes necessary for glucose assimilation into the TCA cycle (KEGG pathway map ID bteq00020) by glycolysis (bteq00010). Likewise, genes for the metabolism sucrose and starch (bteq00500), and fructose (via fructose-1P and frutokinase (bteq00051)) were present. Pathways for galactose, maltose, and raffinose catabolism were also predicted (bteq00052, bteq00500). However, not all the genes needed for inositol utilization were present (bteq00562), such as those that encode transporters (iolF, iolT). Genes needed for scyllo-inosose transformation (iolE, iolD, iolB, iolC, iolJ) were also absent (Supplementary Table S6).

In addition to sugars, organic acids (e.g., l-lactic acid, citric acid, and l-malic acid) are common on the foliage of plants (Morgan and Tukey, 1964). EA-CB0015 utilized these and other compounds as carbon sources in API 50 CHB/E (Villegas-Escobar et al., 2013) and Biolog GEN III Microplate growth assays (Supplementary Table S7). Although methanol is readily available on plant surfaces (Kutschera, 2007), the absence of genes for methanol dehydrogenase suggest EA-CB0015 is unable to consume it as carbon source. However, genes for formaldehyde fixation (hxlA, hxlB) shared high sequence identity (>90%) with their corresponding homologs in 168, suggesting formaldehyde may be assimilated through the ribulose monophosphate pathway (bteq_M00345, Supplementary Table S6). This may allow EA-CB0015 to scavenge formaldehyde produced by neighboring methylotrophs during growth on plant surfaces.

EA-CB0015 contained genes for β-glucanases (bglC, bglS), chitinase (Csn), extracellular proteases (aprE, nprE), bacillopeptidase F protein degradation (Bpr), endolevanase (levC), xylanase (xynA, B, and C), α-amylase (amyE, amyX), and pectate lyases (pel, pelB; Supplementary Table S8). These enzymes are be used to liberate sugars from complex polysaccharides that are then be assimilated by the strain. Chitinase may also function to inhibit nearby fungi by hydrolyzing their cell walls and inhibiting formation of hyphae (Kobayashi and Ikemoto, 2019; Legein et al., 2020). Lastly, a gene encoding for the quorum quenching enzyme YtnP (lactonase-homolog protein) was also found in the genome. This may degrade quorum sensing autoinducers of competitive strains in the same plant environment (Schneider et al., 2012).

EA-CB0015 contained genes required to reduce nitrate to ammonia via the dissimilatory nitrate reduction pathway (bteq_M00530, Supplementary Table S7). The strain lacks genes for a nitrogenase enzyme complex, and thus is unable to fix atmospheric dinitrogen. Genes required for metabolism of aspartate/glutamate (bteq00250), arginine/proline (bteq00330), glycine/threonine/serine (bteq00260) were present (Supplementary Table S7), all of which have been detected in the foliage of plants (Morgan and Tukey, 1964; Parangan-Smith and Lindow, 2013).

Under phosphate limiting conditions, many microorganisms can solubilize inorganic sources of phosphate or mineralize organophosphorous compounds (Hanif et al., 2015). Genes encoding for the pho regulon, including alkaline phosphatase, extracellular enzymes to catabolize organophosphorous compounds (PhoA, PhoB, PhoD, GlpQ), phosphate transporters, and the tat secretion system for protein export were indeed present within the genome (Supplementary Table S9; Allenby et al., 2005). Interestingly, phytase (phy) was only present as a pseudongene (Supplementary Table S9), suggesting strain is unable to hydrolyze phytate, a major source of organic phosphate.

EA-CB0015 encoded genes for the biosynthesis of the siderophore bacillibactin (dhbACEBF), transporters for importing it into the cell (FeuABC-YusV), and the esterase BesA to release the iron in the cytosol (Supplementary Table S10; Miethke et al., 2006). Genes encoding hydoxamate (fhuBCDG), petrobactin/catecholate (fpbOPQ), schizokinen/anthrobactin (yfhAYZ) transporters suggest the strain may scavenge siderophores produced by other species (Supplementary Table S10). Genes for citrate-iron transporter (fecCDEF), high-affinity iron transporter (efeUOB), and heme degradation (hmoAB) suggest these as additional mechanisms for iron acquisition.

Residing on plant leaves and surfaces, epiphytes are continually challenged by photooxidative stress. In addition to sporulation, biofilm formation, and motility, EA-CB0015 contained genes to protect itself from desiccation, UV light, and oxidative damage (Mols and Abee, 2011). The flavin-dependent photoreceptor (ytvA) activates general stress response mechanisms in the presence of blue light. This includes spore-product photolyase (splB; Supplementary Table S11) to repair thymine dimer adducts produced from UV radiation (Herrou and Crosson, 2011; Vanhaelewyn et al., 2020). Genes for catalase, superoxide dismutase, thioredoxin reductase, hydroperoxide reductase, peroxiredoxin, and antioxidants including bacillithiol, phytoene, and sporulene (Supplementary Table S11), may help EA-CB0015 reduce damage from oxidative stress (Engelmann and Hecker, 1996; Newton et al., 2009; Zuber, 2009; Jeong et al., 2018). Production of extracellular polysaccharides and the siderophore bacillibactin (a UV-B absorbing compound) may further allow EA-CB0015 to avoid damage from UV exposure (Supplementary Table S11).

Epiphytic bacteria scavenge nutrients, influence development of plants, and the composition of local microbial communities through their natural products (Legein et al., 2020). To identify the repertoire of molecules EA-CB0015 may produce, we first analyzed its genome using antiSMASH. We found seven BGCs for known natural products and two of unknown function (Supplementary Table S12). Encoded within the genome were non-ribosomal peptide synthetase (NRPS) BGCs for lipopeptides (plipastatin/fengycin C, surfactin, iturin) and the siderophore bacillibactin. The BGCs for subtilosin A (thiopeptide), sporulation killing factor (SKF; sactipeptide), and bacilysin were also detected (Figure 4). Annotation of the NRPS domains for fengycin C (Villegas-Escobar et al., 2013), suggests the strain actually produces plipastatin. Specifically, the lack of an epimerization domain in module 3 of fenB suggests incorporation of l-Tyr instead of d-Tyr. Moreover, the domain annotation of module 9 predicts epimerization, suggesting d-Thr instead of a l-Thr in the final peptide natural product (Hussein, 2019). The neighborhoods from EA-CB0015 were all conserved to reference BGCs for these natural products (Supplementary Figure S4).

Two genomic regions were annotated to encode for potential terpene and polyketide compounds. The first region (2,107,965 to 2,129,124 bp) contained a gene annotated as squalene-hopene cyclase, but its neighborhood lacked commonly associated genes for condensation of farnesyl diphosphate to squalene (SQase, or hpnD, hpnC, hpnE; Supplementary Table S13; van der Donk, 2015). The second region (2,174,022 and 2,213,530 bp) contained a putative type III polyketide synthase (chalcone synthase; BpsA) and an isoprenylcysteine carboxyl methyltraferase (BpsB), suggesting a role for biosynthesis of aliphatic polyketides, such as triketide pyrones, tetraketide pyrones and alkylresorcinols (Supplementary Table S14; Nakano et al., 2009).

Putative pathways for other natural products were identified by KEGG and manual annotation. EA-CB0015 encoded genes for the antibiotic kanosamine, which is biosynthesized in three steps from glucose-6-phosphate (enzymes NtdC, NtdA, NtdB; KEGG accession bteq00998, Supplementary Table S8). Biosynthesis of the plant hormone IAA included genes for tryptophan aminotransferase (patB), idole-3-pyruvate decarboxylase (yclC), and indole-3-acetaldehyde dehydrogenase (dhaS). These would produce indole 3-pyruvic acid pathway (IPyA) the predominant precursor to IAA (Shao et al., 2015). Genes related to the tryptamine-TAM pathway (bsdC, flavin monamine oxidase, and dhaS), known as an alternative pathway for the synthesis of IAA, were found (Shao et al., 2015). YhcX, predicted to act as a nitrilase in the last step of the indole 3-acetonitrile (IAN) pathway, was also present (Supplementary Table S15).

Leaf surfaces are proposed as hot spots for lateral gene transfer and important breeding grounds for microbial diversity (Lindow and Leveau, 2002). Thus, we sought to understand both the degree to which EA-CB0015 may have been affected by gene transfer events and mechanisms it may have to maintain genome integrity.

EA-CB0015 encoded genes needed to take up DNA by genetic competence (Supplementary Table S5), but did not contain plasmids nor any complete integrative or conjugative transposons (Auchtung et al., 2016). However, some pseudogenes (immA, immR, phrI) or integrases (int) for excision of ICEBs were present in the genome (Supplementary Table S16). Four prophage regions were predicted, suggesting EA-CB0015 was susceptible to bacteriophages (Figure 5; Supplementary Figure S5). PHASTER analysis classified three as incomplete (Bacillus phi4J1, Brevibacillus Jimmer2, and Staphylococcus SPβ-like) and one as questionable (Bacillus SPβ). While insertions and recombination from prophages may disrupt important genes, they may also introduce phage resistance or prototrophy to improve the competitive fitness (Kohm and Hertel, 2021).

Other hallmarks of horizonal gene transfer include GEIs and transposable elements. GEIs encode cluster of genes for specialized functions (pathogenicity, symbiosis, metabolism, resistance, fitness) presumed to be of foreign origin (Juhas et al., 2009). Interestingly, GEIs accounted for 10% of the total genome, with many phage-related genes overlapping the same regions (Figure 5; Supplementary Figures S5, S6). Fourteen ISs (transposable elements) were present in the genome (Supplementary Table S17). However, ISs belonging to IS1182 and IS1595 families were pseudogenes and only those belonging to the IS3 family were complete. The latter was the most abundant, with 8 copies of ISBsuI and 1 copy of ISErh (Siguier et al., 2006; Figure 5). These may offer a selective advantage by accelerating genome rearrangement, introducing virulence factors, or resistance genes to antibiotics (Lindow and Leveau, 2002). Genes for putative virulence factors including hemolysin III (hlyIII), capsular polyglutamate and bacillibactin were encoded outside of mobile elements (Figure 5; Supplementary Table S18). Resistance genes for streptomycin (aadK), spiramycin and telithromycin (mphK), and tetracycline (tetL; Figure 5; Supplementary Tables S18, S19) were also detected. The streptomycin resistance gene was encoded within a genomic island.

Both type I and type II RM systems were encoded in the genome of EA-CB0015. The type I system contained genes for the HsdR endonuclease (R), specificity subunit (S), and DNA methyltransferase (M; Figure 5; Supplementary Figure S7). While genes for the endonuclease and methyltransferase of the type I system were highly conserved (>85% sequence identity), the specificity subunit (G4P54_RS16760) exhibited significantly lower sequence identity (33%–40%; Supplementary Figure S8). This suggests variability in the recognition site where the restriction endonuclease cleaves DNA. Additionally, a type II RM system with a site-specific DNA methyltransferase (locus tag G4P54_RS15065) was found.

EA-CB0015 contained 22 TA systems (Figure 5; Supplementary Figure S9; Supplementary Table S20). These included two type I TA systems (bsrG/SR4 and yonT-yoyJ/SR6 located in SPβ region) and three type II TA systems (ndoA/ndoB, spoIISA/spoIISB/spoIISC and yeeF/yezG) which are also in strain 168. TA systems may provide fitness advantages when the strain experiences stress from pH changes, oxygen deficiency, or iron limitation (Brantl and Müller, 2019). Lastly, no RNA-based defense systems were found using CRISPRCasFinder (Couvin et al., 2018).

To understand how EA-CB0015 compares with these closely related bacteria and the presence of any distinguishing trends, we cataloged natural product biosynthesis, prophage, and RM across type strains within the B. subtilis species complex, B. cereus species complex and B. coagulans (Figure 2). In general, strains that are non-pathogenic to humans (B. subtilis, B. tequilensis, B. mojavensis, B. atrophaeus, B. velezensis, B. amyloliquefaciens, B. licheniformis and B. pumilus) encoded a more natural product BGCs. Among the 22 different natural product BGCs we identified, the largest number (7–9) were associated with B. subtilis, B. velezensis, B. amyloliquefacies, and B. tequilensis (Figure 2B). These was consistent with known biological control agents produced by these species (Chen et al., 2009; Gutierrez-Monsalve et al., 2015; Pandin et al., 2018; Samaras et al., 2021). Surfactin, bacilysin, bacillibactin, and bacillaene were the most prevalent BGCs. Except for the antifungal phosphonopeptide rhizocticin (Borisova et al., 2010), EA-CB0015 and BAA-819 share the same composition of natural product BGCs. The most common BGC in pathogenic bacilli was petrobactin. None of these 22 BGCs were present in B. coagulans.

Prophages account for substantial genetic variation and confer phage resistance (Casjens, 2003; Fortier and Sekulovic, 2013; Kohm and Hertel, 2021). Across all taxa, the number of candidate prophage regions and prophage CDSs was similar within the B. subtilis species complex (4.68 regions and 217.74 CDSs) and B. cereus group (5.15 regions and 196.95 CDSs; Figure 2C; Supplementary Figure S10). Intact prophages included members of the Siphoviridae, Myoviridae and Tectiviridae families (Supplementary Figure S11). Twelve phages were unique to B. subtilis species complex, 8 unique to the B. cereus species complex, and 3 were shared between both groups (Supplementary Figures S10, S12). Altogether, these results indicate phages are evolutionary conserved within species complexes.

As natural selection is one of the most important evolutionary processes, RM systems serve as an important defense against the introduction of foreign DNA (Bourniquel and Bickle, 2002). Type II RM systems were the most common among both groups, present in 68 and 81% of the strains from B. subtilis and B. cereus complexes, respectively. Other types of RM systems were also widely distributed but less abundant than type II systems (Figure 2D), but no differential patterns were observed between these two groups of strains.