Characterization and complete genome analysis of the surfactin-producing, plant-protecting bacterium Bacillus velezensis 9D-6

Isolation of bacteria with antimicrobial activities

Out of six hundred bacterial isolates from a remediated potato rhizosphere in Norfolk County, Ontario, Canada, isolate 9D-6 exhibited the highest degree of microbial growth inhibition. Zones of microbial inhibition on solid medium around discs inoculated with 9D-6 indicated that it inhibits growth of the Gram positive bacteria Bacillus cereus and Clavibacter michiganensis; and of the Gram negative bacteria Pantoea agglomerans, Ralstonia solanacearum, Xanthomonas campestris, and Xanthomonas euvesicatoria; but not Erwinia amylovora or Pseudomonas syringae DC3000. It inhibited growth of the ascomycete fungi Cochliobolus carbonum, Fusarium oxysporum, Fusarium solani, Gibberella pulicaris, Gibberella zeae, Monilinia fructicola, and Pyrenochaeta terrestris; the basidiomycete fungus Rhizoctonia solani; and the deuteromycete fungus Alternaria solani; but not the oomycete Pythium mamillatum. In summary, isolate 9D-6 inhibited growth of six of the seven bacterial pathogens and all nine of the fungal pathogens that were tested (Figs. 1 and 2).

We further tested growth of strain 9D-6 at various temperatures and pH. Strain 9D-6 grew in liquid LB medium at a temperature range of at least between 12 °C and 50 °C, and at a pH between 5 and 8, with optimal growth around 30 °C and pH 7 (Additional files 1 and 2).

B. velezensis 9D-6 phylogeny

BLASTn analysis of 16S rDNA amplified from isolate 9D-6 returned a > 99% identity match with the reference strain B. velezensis G341. Therefore, isolate 9D-6 is named B. velezensis 9D-6.

To estimate the evolutionary position of B. velezensis 9D-6 relative to 27 other B. velezensis strains for which complete genomes are available, a phylogenetic tree was constructed based on the RNA polymerase β subunit (rpoB) gene. Compared with the 16S rRNA gene, rpoB provides better resolution between closely related organisms [9]. The phylogenetic analysis suggests that B. velezensis 9D-6 is most closely related to B. velezensis AS43.3, and is also closely related to B. velezensis strains CC09, SB1216, TrigoCor1448, UCMB5033, and UCMB5036, followed together by UCMB5113 and the reference strain G341 (Fig. 3).

In silico DNA-DNA hybridization (isDDH) was performed to confirm 9D-6 as a member of B. velezensis and to validate its position in the phylogenetic tree. Of the strains included in the tree (Fig. 3), the highest isDDH values with the genome of B. velezensis 9D-6 were obtained for strains TrigoCor1448 (92.5%), AS43.3 (92.0%), UCMB5113 (91.7%), G341 (91.7%), UCMB5033 (91.2%), SB1216 (91.1%), UCMB5036 (89.3%), and CC09 (89.1%). These eight strains form a cluster with B. velezensis 9D-6 in the rpoB phylogenetic tree. The values for isDDH between B. velezensis 9D-6 and that of the other B. velezensis strains range from 79 to 86%. In contrast, the value for isDDH is 21% between B. velezensis 9D-6 and B. subtilis subtilis 168. However, it ranges between 55 and 92% between B. velezensis 9D-6 and 22 B. amyloliquefaciens strains with complete genomes: 15 of these values are between 79 and 92%, while the remaining six (including the type strain DSM 7) are between 55 and 56%.

B. velezensis 9D-6 produces antimicrobial compounds [Leu⁷] surfactin C₁₄ and [Leu⁷] surfactin C₁₅

To identify the antimicrobial compounds that contribute to biocontrol, B. velezensis 9D-6 metabolites were isolated into fractions using liquid chromatography, and each fraction was tested for antibacterial activity against the Gram-positive bacterium C. michiganensis, chosen due to its rapid growth rate and agricultural relevance [10, 11]. A fraction demonstrating activity (Fig. 4a) was then subjected to LC-MS to identify the active compound.

The two major peaks observed had predicted chemical formulas of C₅₂H₉₁N₇O₁₃, and C₅₃H₉₃N₇O₁₃. Antibase identified the lipopeptide surfactins as possible identifications. Comparison of the MS/MS of compound C₅₂H₉₁N₇O₁₃ with previously published MS/MS spectra [12] identified this peak as [Leu⁷] surfactin C₁₄ (sometimes called surfactin B). MS/MS of the peak C₅₃H₉₃N₇O₁₃ allowed partial amino acid sequence to be determined by de novo sequencing (Fig. 5), suggesting it is [Leu⁷] surfactin C₁₅ (sometimes called surfactin C).

The surfactin fraction was then tested against the Gram-negative bacterium P. syringae DC3000 and the ascomycete fungus C. carbonum for verification (Fig. 4b-c). The antimicrobial disc inhibited C. carbonum when inoculated with surfactin concentrations between 1 mg/mL and 10 mg/mL. However, for P. syringae DC3000, the zones of inhibition were only apparent at the higher concentrations (5 mg/mL and 10 mg/mL) and were not very large, consistent with the inability of live B. velezensis to antagonize this bacterium (Fig. 1).

B. velezensis 9D-6 reduces the attachment and infection of Arabidopsis by P. syringae DC3000

To determine whether B. velezensis 9D-6 is able to inhibit microbial phytopathogenesis by mechanisms other than antimicrobial production, we tested its ability to reduce root colonization by P. syringae DC3000, a Gram-negative bacterium. Since growth of P. syringae DC3000 was not inhibited by B. velezensis 9D-6 in the in vitro plate assay (Fig. 1), reduced root colonization in response to B. velezensis 9D-6 should be due to a mechanism other than production of bacteriocidal or bacteriostatic antimicrobials.

Co-inoculation of A. thaliana seedlings with B. velezensis 9D-6 and P. syringae DC3000 produced milder disease symptoms compared to inoculation with P. syringae DC3000 alone. By contrast, B. velezensis 9D-6 alone produced no visible symptoms of disease (Fig. 6). Co-inoculation of A. thaliana seedlings with B. velezensis 9D-6 and RFP-labeled P. syringae DC3000 also resulted in a notable decrease in root colonization by the pathogen at all tested concentrations, as indicated by an observable decrease in RFP fluorescence (Fig. 7).

The B. velezensis 9D-6 genome sequencing and analysis

To allow further investigation into the biocontrol mechanisms of B. velezensis 9D-6 and its application for sustainable agriculture, we determined its complete genome.

The B. velezensis 9D-6 genome consists of a single 3.96 Mb circular chromosome (Fig. 8), which fits within the range of 3.81 to 4.24 Mb reported for other completed genomes of the species. Like nearly 80% of those strains, it does not harbor a plasmid. The B. velezensis 9D-6 genome is predicted to include 3942 total genes, of which 3849 (97.6%) are protein coding genes, 93 (2.4%) are RNA coding, and 82 (2.1%) are pseudogenes (Additional file 3). Of the RNA coding genes, 21 are predicted to code rRNA and 68 are predicted to encode tRNA. Among the predicted genes, 2736 (69.4%) are associated with general COG function categories. The distribution of genes in these categories is presented in Additional file 4.

A core genome containing 2574 coding sequences is shared among B. velezensis, the closely related species Bacillus amyloliquefaciens, and B. subtilis subtilis 168 [13]. Among B. velezensis strains with completed genomes, the number of coding sequences is typically between 3600 and 4000 (with a clear outlier of 4411 for B. velezensis SRCM101413) and the number for B. velezensis 9D-6 fits within this expected range. Similarly, the 46.4% GC content of B. velezensis 9D-6 is comparable to the other completed B. velezensis genomes, which range from 45.9 to 46.8% GC. The 9D-6 genome sequence was deposited in GenBank (CP020805): https://www.ncbi.nlm.nih.gov/nuccore/CP020805

Special genomic features

Manual mining of B. velezensis 9D-6 also identified a gene for butanediol dehydrogenase and genes predicted to encode components of secretion systems including the inner-membrane spanning twin-arginine translocation (Tat) system, SecYEG translocon, and the Type VII/ESX secretion system.

Isolation of B. velezensis 9D-6

Soil samples were collected from Blizman potato fields in Norfolk County, Ontario, Canada in the summer of 2012. Over the previous three years, bio-organic fertilizer was added to the soil each spring in an effort toward natural remediation. In 2012 (the fourth year), 10.0 g of moist soil was collected, placed in 95 mL of sterile water, and shaken for 10 min. Then, 1.0 mL of this suspension was transferred for serial dilution up to 10^− 10, and the dilutions were plated on tryptic soy agar (TSA) for 48 h at 28 °C to attain single colonies. Permission for this research was obtained from the Canadian Food Inspection Agency (CFIA).

In vitro antagonism assays against microbial pathogens

Six-hundred bacterial strains were further screened for their ability to suppress phytopathogens under in vitro conditions. Antimicrobial discs were prepared by inoculating 0.5 mm discs of P8 Filter Paper (Thermo Fisher Scientific, Pittsburgh, PA, USA) with 50 μL of 10⁹ CFU/mL B. velezensis 9D-6 suspended in 0.85% NdaCl. To test inhibition of pathogenic bacteria, 100 μL of each bacterial test strain at 10⁹ CFU/mL in 0.85% NaCl were spread onto separate LB agar plates, and an antimicrobial disc containing B. velezensis 9D-6 was placed in the centre of the plate. After two days incubation at 28 °C, zones of inhibition around the discs were recorded. Bacterial test strains were Bacillus cereus, Clavibacter michiganensis, Erwinia amylovora, Pantoea agglomerans, Pseudomonas syringae DC3000, Ralstonia solanacearum, Xanthomonas campestris, and Xanthomonas euvesicatoria.

To test inhibition of pathogenic eukaryotes, four plugs from each tested strain were evenly spaced at the peripheries of a potato dextrose agar plate, and an antimicrobial disc containing B. velezensis 9D-6 was placed in the centre of the plate. After one to two weeks of incubation (depending on the growth rate of the test strain) at room temperature, zones of inhibition around the discs were recorded. Eukaryotic test strains were Alternaria solani, Cochliobolus carbonum, Fusarium oxysporum, Fusarium solani, Gibberella pulicaris, Gibberella zeae, Monilinia fructicola, Pyrenochaeta terrestris, Pythium mamillatum, and Rhizoctonia solani.

16S rRNA sequencing and analysis

Bacterial genomic DNA was used for PCR amplification of 16S rRNA. For genomic DNA isolation, a single bacterial colony was inoculated into 2.5 ml of LB broth and grown for 16 h at 30 °C with shaking at 200 rpm. Cells were collected by centrifugation of 1.5 ml culture at 13,000 rpm for 5 min and bacterial genomic DNA was isolated using GenElute Bacterial Genomic DNA Isolation kit (Sigma-Aldrich Co., St. Louis, MO, USA) in accordance with the manufacturer’s protocol. PCR amplification of an approximately 1,500 base pair sequence of the bacterial 16S rDNA gene was performed with primers 8F (5′- AGAGTTTGATCCTGGCTCAG-3′) and 1492R (5’-GGTTACCTTGTTACGACTT-3′) [35]. Each 50 μL PCR mixture contained 1.5 units Phusion High-Fidelity DNA Polymerase (Thermo Fisher Scientific Inc., Waltham, MA, USA), 1X PCR buffer, 2.0 mM MgCl2, 200 μM dNTPs, 2.5 μM of each primer, 50 ng of genomic DNA template, and ultrapure water (Sigma-Aldrich Co.). The cycle parameters were as follows: initial denaturation at 95 °C for 5 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 57 °C for 45 s, and extension at 72 °C for 60 s, with the final overall extension at 72 °C for 10 min. The 16S PCR products were purified by using QIAquick PCR Purification kit (Qiagen, Hilden, Germany) and sequenced with 8F and 1492R primers on a 3730 DNA Analyzer (Thermo Fisher Scientific Inc.) at Agriculture and Agri-Food Canada (London, Ontario, Canada). The 16S rRNA gene fragment was compared with the NCBI nucleotide database using Blastn to determine the closest taxonomic relatives.

Phylogenetic tree construction

A phylogenetic tree was constructed based on RNA polymerase β subunit (rpo) gene sequences [9] obtained from the genome sequence (see below) to estimate the evolutionary position of B. velezensis 9D-6 relative to other B. velezensis strains. The Rpo protein sequence of Bacillus subtilis subspecies subtilis strain 168 (GenBank Accession NP_387988) was used as a query for translated BLAST (tblastn) against B. velezensis gene sequences in the NCBI database. In total, 28 B. velezensis rpo sequences were copied and aligned with rpo from B. subtilis subtilis 168 using WebPRANK software [36]. The type strains of B. velezensis are NRRL B-41580 and KCTC 13012 [37], but neither of their genome sequences are complete, and they were therefore excluded from the analysis. Phylogenetic trees were constructed using MrBayes v.3.2 [38] with the generalised time reversible (GTR) substitution model and default parameters for run length and frequencies. The B. subtilis subtilis 168 rpo was designated as an outgroup. Following analysis, the standard deviation of spilt frequencies was 0.006566 and the potential scale reduction factor (PSRF) was 1.000 or 1.001 for all factors. The majority-rule consensus tree was visualized using TreeGraph2 [39].

Extraction and identification of antimicrobial compounds

To identify potential antimicrobial compounds produced by B. velezensis 9D-6, metabolites were isolated from a liquid LB culture grown at 28 °C for 72 h by ethyl acetate, liquid/liquid extraction. The extract was dried and resuspended in acetonitrile. After verifying the antibacterial and antifungal activity of the crude solution, compounds were separated into 22 fractions using a C18 reverse-phase high performance liquid chromatographer on a 1260 Infinity Series (Agilent Technologies, Santa Clara, CA, USA) eluted using an acetonitrile-water gradient with 0.1% formic acid.

Fractions were each tested in duplicate for their antibacterial activity against C. michiganensis. Each antimicrobial disc was prepared by adding 50 μL of a fraction to a 0.5 mm disc of P8 Filter Paper (Thermo Fisher Scientific, Pittsburgh, PA, USA), which was then placed onto an LB agar plate spread with 100 μL of 109 CFU/mL C. michiganensis in 0.85% NaCl. Controls used acetonitrile in place of the metabolite fraction.

Fractions demonstrating antibacterial effects were individually analyzed using ultraviolet-visible (UV-vis) spectroscopy at 210 nm and 254 nm to identify UV-active compounds. These compounds were isolated and tested in a final in vitro antibacterial plate assay against C. michiganensis. The compounds of bioactive fractions were characterized by high resolution LC-MS/MS using a Q-Exactive Quadrupole-Orbitrap Mass Spectrometry (Thermo Fisher Scientific, Waltham, MA, USA). Full MS and MS/MS spectra were analyzed manually with XCalibur software (Thermo Fisher Scientific, Waltham, MA, USA). The full MS data was used to identify possible molecular formula (< 3 ppm). Formulas were than searched against the AntiBase 2012 Natural Compound Identifier database [40] to identify purified compounds. When possible, the MS/MS spectra were used for de novo sequencing of peptide containing compounds to confirm possible Antibase identified compounds.

The surfactin fraction was then tested for activity against P. syringae DC3000 and C. carbonum. Discs were inoculated with 50 μL of the fraction at concentrations of 10 mg/mL, 5 mg/mL, 2.5 mg/mL, and 1 mg/mL. Discs were then placed onto LB agar plates that were either spread with P. syringae DC3000 or contained a central agar plug of C. carbonum.

In vivo biocontrol

To determine whether alternative control mechanisms (such as induced systemic resistance) might be elicited by B. velezensis 9D-6, we tested its ability to reduce root colonization by P. syringae DC3000, a Gram-negative bacterium whose growth in vitro was not inhibited by B. velezensis 9D-6. P. syringae DC3000 was labeled using a red fluorescent protein (RFP) reporter gene construct based on the plasmid pME6010 [41, 42], and used to inoculate Arabidopsis thaliana (L.) Heynh seedlings in hydroponic medium (2.5 × 10⁵ CFU/mL), with or without B. velezensis 9D-6 (2.5 × 10⁶ CFU/mL), as previously described [43]. The A. thaliana (L.) Heynh seeds were obtained from the Arabidopsis Biological Resource Center (Columbus, OH, U.S.A) and grown as previously described [43].

After seven days, A. thaliana roots were removed from the medium, rinsed in ultrapure water to remove loosely bound material, and attachment by labeled P. syringae DC3000 was imaged using a DMIRE2 inverted microscope with confocal laser scanner (Leica Microsystems GmbH, Wetzlar, Germany). Samples were excited using a helium-neon 543/594 nm laser, and emission was detected at 590–630 nm under a 63 x water immersion objective with a numerical aperture of 1.4.

Genome sequencing and annotation

Following isolation from potato field soil of Norfolk County, Ontario, Canada, a single colony of B. velezensis 9D-6 was grown for 48 h in 3 mL of lysogeny broth (LB) at 28 °C and 60 rpm in a TC-7 drum rotor (New Brunswick Scientific Co., Enfield, CT, USA). Genomic DNA was isolated from 1.5 mL of the culture using the GenElute Bacterial Genomic DNA kit (Sigma-Aldrich Co., St. Louis, MO, USA) according to the manufacturer’s protocol for Gram positive bacteria, with the exception that DNA was eluted in UltraPure DNase/RNase-free distilled water (Thermo Fisher Scientific Corp., Waltham, MA, USA). DNA quality was assessed by agarose gel electrophoresis.

Genomic DNA was sequenced at AGCT Inc. (Wheeling, IL, USA) using the MiSeq next generation sequencing platform (Illumina Inc., San Diego, CA, USA). Libraries were constructed using the NexteraXT DNA sample preparation kit (Illumina Inc., San Diego, CA, USA) with a target average insert size of 500–600 bp. Sequencing generated 5,453,964 raw read pairs of 2 × 150 bp read length on average. Adaptor sequences, low quality sequences, and short reads were filtered out using BaseSpace (Illumina Inc., San Diego, CA, USA), Trim Galore (Babraham Bioinformatics, Babraham, UK), and Sickle [44]. The remaining 4,722,973 trimmed read pairs were assembled de novo and scaffolded with SPAdes v. 3.50 [45] into 47 contigs with an N₅₀ of 507 kb spanning 100x coverage.

The draft genome was aligned with test reference genomes, including that of B. velezensis G341 (GenBank Accession CP011686), using the Mauve Multiple Genome Alignment tool [46] version 2.4.0 to generate a preliminary chromosome map. Primers were designed based on the contig orientations suggested by the chromosomal map, and connections were confirmed by Sanger sequencing of PCR products followed by alignment with the draft genome using SeqMan Pro (DNASTAR Inc., Madison, WI, USA) [47]. All gaps were closed with one of two repeat sequences; one containing a transposon region and one containing an rDNA operon sequence. After closure and validation of gaps between draft contigs, the final completed genome was assembled into a single 3.96 Mb chromosome using the largest 20 contigs. The remaining 27 contigs were identified as redundant fragments.

Numbers of genes (rRNA genes, protein coding genes, pseudogenes) were predicted on the DOE-JGI Microbial Genome Annotation Pipeline (MGAP v.4) [48]. In silico DNA-DNA hybridization (isDDH) was performed via the GGDC web server (http://ggdc.gbdp.org) using formula 2 [49]. Secondary metabolite clusters were identified using both MGAP v.4 and antiSMASH3.0 [50]. The genome was also mined in silico using target amino acid sequences from closely related species from the UniProtKB sequence database [51] and comparing them to translated draft genome nucleotide data using tBLASTn.