Draft genome sequence of Bacillus velezensis 2A-2B strain: a rhizospheric inhabitant of Sporobolus airoides (Torr.) Torr., with antifungal activity against root rot causing phytopathogens.

A Bacillus velezensis strain from the rhizosphere of Sporobolus airoides (Torr.) Torr., a grass in central-north México, was isolated during a biocontrol of phytopathogens scrutiny study. The 2A-2B strain exhibited at least 60% of growth inhibition of virulent isolates of phytopathogens causing root rot. These phytopathogens include Phytophthora capsici, Fusarium solani, Fusarium oxysporum and Rhizoctonia solani. Furthermore, the 2A-2B strain is an indolacetic acid producer, and a plant inducer of PR1, which is an induced systemic resistance related gene in chili pepper plantlets. Whole genome sequencing was performed to generate a draft genome assembly of 3.953 MB with 46.36% of GC content, and a N50 of 294,737. The genome contains 3713 protein coding genes and 89 RNA genes. Moreover, comparative genome analysis revealed that the 2A-2B strain had the greatest identity (98.4%) with Bacillus velezensis.

Introduction

Root rot causing microorganisms are among the most devastating phytopathogens of many horticultural crops resulting in considerable financial loss worldwide. Some of these pathogens include the oomycete , and the fungi , , and . Biocontrol strategies are important alternatives to keep some plant pathogens at low levels in affected crops, particularly when evaluating the risk of the use of pesticides on human health and the environment, and the social pressure to have innocuous horticultural food products. Biocontrol agents including bacterial strains that possess biocide activity against phytopathogens can also have the ability to invoke a systemic resistance (induced systemic resistance) in the host plant [1, 2]. In some cases, these bacterial strains are also able to promote the plant growth by inducing the biosynthesis of phytohormones [3, 4]. The rhizosphere is the area around the plant root that is inhabited by a unique population of microorganisms. The rhizospheric space is characterized by plant root exudates and usually by a high density and diversity of microorganisms. The root exudates have positive and negative effects in the interactions in the rhizosphere [5-7]. Colonizers of the rhizosphere are a great variety of microorganisms including bacteria that commonly have a friendly interaction with the plant host, suppressing at the same time some phytopathogens, and in some cases promoting plant growth [2, 7].

In seeking new options for biocontrol alternatives against phytopathogens, the genomic and biotechnological advances allows the deciphering of the molecular processes that regulate and induce the expression of many genes of plant-associated microorganisms. This will increase the possibilities for newer options of biocontrol agents with improved efficacy to deal with specific pathological problems of important crops.

In the present study we sampled soil and roots of cultivated and wild plants from Zacatecas state in the central-north region of Mexico, and isolated the bacterial strains from rhizosphere of (Torr.) Torr. The bacteria reported in this study has the capacity to inhibit the growth of each of the four virulent isolates of the pathogens , , and , which are the causal agents of root rot in chili pepper crops. Levels of indoleacetic acid and biosynthesis of siderophores were analyzed, along with the induction of NPR1, a key gene controlling local resistance and systemic acquired resistance with multiple roles in plant immunity [8, 9]. In addition, the expression of the sesquiterpene cyclase gene involved in the isoprenoids pathway for the biosynthesis of phytoalexin capsidiol [10] was analyzed. The bacterial strain 2A-2B from this rhizosphere was selected for genome sequencing. A draft genome assembly of 3.953 MB was obtained and deposited in the NCBI GenBank (biosample SAMN05772828 and accession MLCV00000000), and in the Genomes OnLine Database (GOLD) with the accession Gp0177877.

Organism information

Classification and features

10.1601/nm.8505 strain 2A-2B is a Gram-positive bacterium, with rapid growth rate in LB liquid medium reaching the stationary phase after 13 h at 28 °C. In contrast, the growth rate of this strain was much slower in LB solid medium with the stationary phase attained at 24 h. The colonies in solid LB medium were observed as circular with pulvinated elevation and had filiform beige-opaque margins.

The strain was isolated from the rhizosphere of (Torr.) Torr., a wild grass growing in a grassland area of Morelos municipality in Zacatecas state, Mexico. This bacterium could grow well in mediums with high content of nutrients as LB, KB, TSA, PDA, and also in YMB medium. The optimal temperature for growth of this bacterium is 28 °C, although it was also capable of growing in temperatures of up to 35 °C. In root colonization assays in vitro in “mirasol” pepper plantlets ( L. mirasol), no detrimental effect was seen in the plant. Furthermore, the 2A-2B 10.1601/nm.8505 strain produces indoleacetic acid (data not shown), a plant growth regulating auxin, but does not synthesize siderophores. The morphology of this motile rod shaped Gram-positive bacterium is shown in Fig. 1 and the general features in Table 1.

Fig. 1

Cellular morphology of B. velezensis 2A-2B. Gram-positive bacteria under phase contrast with a Leica DM2500 compound microscope with Nomarski differential interference contrast

Table 1

Classification and general features of Bacillus velezensis strain 2A-2B [37]

MIGS ID	Property	Term	Evidence codea
	Classification	Domain: Bacteria	TAS [38]
		Phylum: Firmicutes	TAS [39, 40]
		Class: Bacilli	TAS [41]
		Order: Bacillales	TAS [42]
		Family: Bacillaceae	TAS [43]
		Genus: Bacillus	TAS [44]
		Species Bacillus velezensis	TAS [37]
		Strain: 2A-2B
	Gram stain	Positive	IDA
	Cell shape	Rod	IDA
	Motility	Motile	IDA
	Sporulation	Not specified	NAS
	Temperature range	20–35 °C	IDA
	Optimum temperature	28 °C	IDA
	pH range; Optimum	6.5–7.0; 6.5	IDA
	Carbon source	Heterotrophic	IDA
MIGS-6	Habitat	Soil	NAS
MIGS-6.3	Salinity	16% (w/v) NaCl	IDA
MIGS-22	Oxygen requirement	Aerobic	IDA
MIGS-15	Biotic relationship	Rhizosphere	NAS
MIGS-14	Pathogenicity	Non-pathogen	IDA
MIGS-4	Geographic location	México/Zacatecas, Morelos.	NAS
MIGS-5	Sample collection	2012	NAS
MIGS-4.1	Latitude	N 22° 49′ 12.354′´	NAS
MIGS-4.2	Longitude	W 102° 41′ 51.59299′´	NAS
MIGS-4.4	Altitude	2200.1 M	NAS

aEvidence codes - IDA: Inferred from Direct Assay; TAS: Traceable Author Statement (i.e., a direct report exists in the literature); NAS: Non-traceable Author Statement (i.e., not directly observed for the living, isolated sample, but based on a generally accepted property for the species, or anecdotal evidence). These evidence codes are from the Gene Ontology project [45]

The phylogenetic analysis based on 16S rRNA sequences using MEGA7 software [11] showed that 10.1601/nm.8505 2A-2B is evolutionarily positioned between 10.1601/nm.8505, 10.1601/nm.4868 and 10.1601/nm.20158 (Fig. 2). In recent studies of genome sequencing and comparative genomics of 10.1601/nm.8505 10.1601/strainfinder?urlappend=%3Fid%3DNRRL+B-41580, 10.1601/nm.20158 10.1601/strainfinder?urlappend=%3Fid%3DKACC+13015 and 10.1601/nm.22414 FZB42, it was established that these last two strains are heterotypic synonyms of 10.1601/nm.8505 [12], in this context and based in our results, the 2A-2B strain corresponds to the 10.1601/nm.8505 species.

Fig. 2

Phylogenetic tree with B. velezensis 2A-2B showing its relationships with species of the Bacillus genus. The tree was constructed based on 16S rRNA gene alignments in MEGA7 [11] applying the Neighbor-join method, and rooted using Bacillus gibsoni FJAT-10019 strain as an outgroup

The 2A-2B strain of 10.1601/nm.8505 inoculated in roots of plantlets of L. mirasol cultivar induced expression of the sesquiterpene cyclase and NPR1 genes (Fig. 3). Sesquiterpene cyclase is involved in the phytoalexin capsidiol biosynthesis pathway in C. annuum [10] and NPR1 gene is a master regulator of systemic acquired resistance in response to biotic stress in plants [8, 13]. This work, in relation to the sampling of materials in the field and the activities performed in laboratory was done under national guidelines.

Fig. 3

NPR1 and Sesquiterpene cyclase genes induced in root inoculated plantlets of Capsicum annuum L. Genes level induction assessed by semi-quantitative RT-PCR. Assays with mixed total RNA from root and leaves of mirasol cultivar pepper plantlets in root inoculated with B. velezensis 2A-2B strain, at 12 hs post-inoculation. Lane 1, control plant not inoculated (mock); lane 2, plant in root inoculated with 2A-2B strain; lane 3, a B. velezensis strain inducer of SC gene in chili pepper. SC, sesquiterpene cyclase; NPR1, nonexpresser of PR genes

Extended feature descriptions

In physiological studies it is of relevance that this 10.1601/nm.8505 2A-2B strain is tolerant to salt. Tests of bacterial growth on solid medium with NaCl show that this bacterium is able to support up to 16% (w/v) (Table 1), whereas in other 10.1601/nm.4857 species and strains, the reported salt tolerance reach up to 12% [14]. This data motivated further studies on plant growth under saline soil condition with 2A-2B strain inoculations. In the genome of this strain five copies of the glycine betaine/L-proline ABC transporter ATP-binding protein are found, and their alignment with Clustal Omega software [15] shows both semiconservative and nonconservative amino acids (Fig. 4). Furthermore, it is found a glycine/betaine ABC transporter permease. This transport system is involved in osmoregulation by accumulating glycine betaine and other solutes under conditions of stress, with the functionality of osmoprotection [16, 17]. The fact that this bacterium contains two genes, one with 5 copies, related to the glycine betaine accumulation, suggests a possible role on the capacity of this microorganism in the osmoprotection mechanism under salt stress conditions.

Fig. 4

Alignment of glycine betaine/L-proline ABC transporter proteins present in the genome of B. velezensis 2A-2B strain. Semiconservative and nonconservative amino acids predominate

Genome sequencing information

Genome project history

The 10.1601/nm.8505 2A-2B strain was selected due to its capacity to inhibit the growth of four pathogens, the causal agents of root rot in chili pepper. These include: , , , and . The 2A-2B strain was not pathogenic in “mirasol” chili pepper when inoculated in the roots.

The genome project of 10.1601/nm.8505 2A-2B was deposited in the GenBank NCBI database as a BioProject, Bio-sample and genome. The IDs for BioProject and genome accession numbers are PRJNA343056 and MLCV00000000, respectively. In the GOLD (genomes on line database) server, the assigned accession number is Gp0177877. A summary of the project information is shown in Table 2.

Table 2

Project information

MIGS ID	Property	Term
MIGS 31	Finishing quality	Quality draft
MIGS-28	Libraries used	illumina paired-end library
MIGS 29	Sequencing platforms	MiSeq illumina
MIGS 31.2	Fold coverage	38
MIGS 30	Assemblers	SPAdes Genome Assembler 3.8.1
MIGS 32	Gene calling method	NCBI Prokaryotic Genome Annotation Pipeline
	Locus Tag	BKK82
	Genbank ID	MLCV00000000
	GenBank Date of Release	2017–02-23
	GOLD ID	Gp0177877
	BIOPROJECT	PRJNA343056
MIGS 13	Source Material Identifier	2A-2B
	Project relevance	Biotechnological, Agricultural, Biocontrol

Growth conditions and genomic DNA preparation

A sample taken from a colony of the 2A-2B strain was inoculated in 5 ml of LB liquid medium and cultured for 24 h at 150 rpm at 28 °C. The bacterial culture was centrifuged at 3500 g and the bacterial pellet was subjected to DNA extraction and purification based on the bacterial DNA isolation CTAB protocol [18]. The purity and concentration of the DNA was analyzed by agarose gel electrophoresis and in a Qubit 2.0 Fluorometer (Invitrogen). One nanogram of 10.1601/strainfinder?urlappend=%3Fid%3DDNA+in+5 ul of water was used to construct the genome libraries tagmented by PCR for Illumina sequencing. The quality and size of fragments in the libraries was verified in a Bioanalyzer (BioAnalizer 2010, Agilent Technologies). The libraries were subjected to standard normalization and 15 pM were used in the sequencing process.

Genome sequencing and assembly

Purified genomic bacterial DNA was used to prepare libraries following Nextera Kit instructions (Illumina, San Diego Ca.). High-throughput sequencing was done under sequencing by synthesis protocol (MiSeq, Illumina) with a 2 × 75 paired-end run at the Sequencing Laboratory at the Unidad de Ciencias Biológicas, Universidad Autónoma de Zacatecas, México. SPAdes Genome Assembler 3.8.1 [19] was used to assemble the genome and the quality of the assembly was evaluated using QUAST 4.1 [20].

Genome annotation

The protein-coding genes, structural RNAs and tRNAs in the draft genome were predicted using the NCBI Prokaryotic Genome Annotation Pipeline [21] and the GENIX Automated Bacterial Genome Annotation Pipeline [22].

Clusters of orthologous groups were predicted using the Web CD-search tool working with the COG database of orthologous protein families focusing on prokaryotes. Genes with Pfam domains and signal peptides were analyzed using the Pfam database [23], and the PrediSi software tool [24] respectively. Genes with transmembrane helices were analyzed with the CRISPRs web server [25]. A summary of all features of genome annotation is shown in Table 3.

Table 3

Genome statistics

Attribute	Value	% of Total
Genome size (bp)	3,958,607	100.0
DNA coding (bp)	3,545,146	89.55
DNA G + C (bp)	1,835,210	46.36
DNA scaffolds	64	100.0
Total genes	3891	100.0
Protein coding genes	3713	95.42
RNA genes	89	2.28
Pseudo genes	89	2.28
Genes in internal clusters	–	–
Genes with function prediction	3388	87.07
Genes assigned to COGs	3026	77,77
Genes with Pfam domains	1329	34.15
Genes with signal peptides	757	19.45
Genes with transmembrane helices	1049	26.95
CRISPR repeats	1	–

Genome properties

The draft genome contains 3,958,607 bp and a GC content of 46.36% with 3713 predicted protein-coding genes. Within this genome, 3388 genes had a predicted function, and 3026 were assigned to clusters of orthologous groups (COGs). In addition to the 3026 genes that were associated to general COG functional categories, 865 genes had not hit in COGs; this information is presented in Table 4, and the map of the genome generated with CGView comparative genomics tool [26] is represented in Fig. 5. Based on the results from sequence annotation and CD-search tool, we found 82 genes in the category of mobilome: prophages and transposons, suggesting that this 10.1601/nm.8505 2A-2B strain could eventually carry extrachromosomal elements when the prophage induction process takes place.

Table 4

Number of genes associated with general COG functional categories

Code	Value	%age	Description
J	221	5.60	Translation, ribosomal structure and biogenesis
A	0	0.0	RNA processing and modification
K	247	6.25	Transcription
L	112	2.84	Replication, recombination and repair
B	0	0.0	Chromatin structure and dynamics
D	63	1.60	Cell cycle control, Cell division, chromosome partitioning
V	78	1.98	Defense mechanisms
T	149	3.77	Signal transduction mechanisms
M	210	5.32	Cell wall/membrane biogenesis
N	61	1.54	Cell motility
U	28	0.71	Intracellular trafficking and secretion
O	112	2.84	Posttranslational modification, protein turnover, chaperones
Z	3	0.07	Cytoskeleton
C	167	4.23	Energy production and conversion
G	250	6.33	Carbohydrate transport and metabolism
E	279	7.07	Amino acid transport and metabolism
F	96	2.43	Nucleotide transport and metabolism
H	169	4.28	Coenzyme transport and metabolism
I	125	3.17	Lipid transport and metabolism
P	145	3.67	Inorganic ion transport and metabolism
Q	99	2.50	Secondary metabolites biosynthesis, transport and catabolism
X	82	2.07	Mobilome: prophages, transposons
R	173	4.38	General function prediction only
S	215	5.44	Function unknown
–	865	21.90	Not in COGs

Fig. 5

Circular map of the chromosome of Bacillus velezensis 2A-2B obtained with CGView comparative genomics tool [26]. From inner to outer rings: ring 1 scale marked in every 200 kbp, ring 2 GC skew (green +, purple -), ring 3 GC% content, ring 4 CDSs on reverse strand and ring 5 CDSs on forward strand

Insights from the genome sequence

Regarding the exhibited antifungal activity of this 10.1601/nm.8505 strain, we found genes in the genome that code for proteins of the Bac operon and the oligopeptide permease OppA. The Bac proteins are involved in the biosynthesis of bacilysin, a non-ribosomally synthesized dipeptide that is active against a range of bacteria and some fungi. The proteolysis of this dipeptide releases the non-proteinogenic amino acid L-anticapsin, which functions as a competitive inhibitor of glucosamine synthase and can result in the lysis of fungal cells [27, 28]. Also, a beta-glucanase and an endoglucanase are present in the genome of this bacterium. Similarly, surfactin synthetase gene which is present in the 2A-2B strain genome, adds to the capacity of this bacterium to contribute in the antifungal activity against the root rot causal agents. Furthermore, the surfactin lipopeptide of 10.1601/nm.10618 is well documented as elicitor of induced systemic resistance in plants [29-31]. In the genome of the 2A-2B strain of 10.1601/nm.8505, with a total of 3713 predicted-protein coding genes, the 1.98% corresponds to defense genes; and the 2.5% of genes corresponds to secondary metabolites biosynthesis. In these two functional categories of genes, a possible role in fungal inhibition may be important. In addition, the sesquiterpene cyclase and NPR1 genes induced in chili pepper plantlets, during the 2A-2B strain root inoculation experiments, suggests that this lipopeptide is sensed by the signaling pathway in the plant’s defense system.

In other hand, in relation to the root bacterial colonization, the CheA and CheY genes are present in the genome of 2A-2B strain. These genes encode proteins that act as a two component system of bacterial chemotaxis, which is a response to chemical signals for controlling the direction of flagellar rotation [32-34]. With this two-component chemotaxis system and other plant exudate chemoreceptors, this bacterium could effectively reach the root tissue and proceeds with the plant tissue colonization.

In the carbohydrate transport and metabolism category, 250 genes (6.3% of total genes) were predicted in the 2A-2B strain genome including the PTS trehalose-specific enzyme IIBC component, xylose isomerase and a number of genes related to glucose metabolism. This suggests that the 2A-2B bacterium possesses a broad battery of genes coding for enzymes required to release a variety of carbon sources including some from plant exudates in the rhizosphere.

Extended insights

The B.velezensis 2A-2B strain has a standard genome size compared to others 10.1601/nm.4857 species where the mean size fluctuate around 3.7 MB; the 2A-2B strain contains a genome of 3.96 MB. It is remarkable that the 2A-2B strain has only one 16S and one 23S rRNAs, whereas other 10.1601/nm.4857 species possesses seven, nine or ten of each.

In genome comparison analysis to another accession of 10.1601/nm.8505 species, the alignment of 2A-2B and 10.1601/nm.8505 FZB42 (acc: NC_009725.1) as reference genome using VISTA software [35], shows blocks of synteny along the reference genome (Fig. 6A); in the genome of 2A-2B strain the great majority of synteny is found in the sense DNA strand, and a minor part in the antisense strand (Fig. 6B). The reference genome shows also fragments of synteny with multiple locations in the compared genome (Fig. 6A).

Fig. 6

Synteny analysis between B. velezensis 2A-2B and B. velezensis FZB42 genomes. The B. velezensis FZB42 (acc: NC_009725.1) as reference genome. In a, synteny of 2A-2B strain vs. FZB42 strain. Blocks of synteny along the reference genome, and remarkable synteny with the sense DNA strand of 2A-2B strain is shown. In b, rearrangements of genomic segments in the 2A-2B strain in comparison to the genome of FZB42 strain. The synteny analysis was performed with VISTA software [35]

In multiple genome comparisons using MAUVE method, which is helpful in the localization of blocks of collinearity, rearrangements and inversions in conserved regions in genomes [36], in the 2A-2B strain of 10.1601/nm.8505 compared to other 6 genomes as a sample from the 89 10.1601/nm.8505 genomes so far accessed in the NCBI Gene Bank, 19 locally collinear blocks (LCBs) were found in the 2A-2B genome, 10 of mayor size, between 28,967 and 808,567 pb length, and 9 of small size, between 208 and 17,724 pb length. Four unique regions in the range of 1010 and 19,421 pb length localized outside of LCBs in the 2A-2B genome, without correspondence in the other 6 10.1601/nm.8505 genomes. Taking as reference any of the other 6 10.1601/nm.8505 genomes, the 6°, 7° and 8° LCBs in the 2A-2B genome are in inverted orientation. The 7° LCB in the 2A-2B genome is deleted in the SCGB574 10.1601/nm.8505 genome, and the 8° LCB of the 2A-2B genome is only conserved among the 2A-2B and 9D-6 of these 7 compared genomes. This alignment of 7 10.1601/nm.8505 genomes shows also that in the 2A-2B genome two events of large rearrangements have occurred, the greater LCB in 2A-2B genome is localized ahead of 4 LCBs in forward orientation compared to the other 6 genomes, and a segment that include three LCBs in the 2A-2B genome is in inverted orientation, in opposed situation in the other 6 compared genomes (Fig. 7).

Fig. 7

Multiple genome comparisons of 2A-2B and other 6 genomes, as a sample of the 89 B. velezensis genomes so far accessed in the NCBI GeneBank. In the 2A-2B B. velezensis strain genome 19 locally collinear blocks (LCBs) are found, this genome with two mayor rearrangements that imply the mayor LCB and a segment including three other LCBs in inverted orientation. The genome comparisons were performed with MAUVE software [36]

Conclusions

In this study, we obtained and characterized a draft genome of 10.1601/nm.8505, the 2A-2B strain isolated from (Torr.) Torr. rhizosphere. The assembled genome contains a total of 3891 genes of which a high number of genes correspond to amino acid and carbohydrate transport and metabolism categories, and transcription and translation functions; whereas the 5.44% of genes have unknown functions. The antagonist characteristics of this bacterium against several fungal phytopathogens, could be explained in part in a number of described genes in the genome that are involved in the biosynthesis of compounds with toxic effects on fungal cell structures.

In comparative analysis with bacterial genomes in the NCBI DataBank using the draft genome and using the16S rRNA sequence, the 2A-2B strain was classified as 10.1601/nm.8505. The genome of the 2A-2B strain contains genes related to antifungal activity and systemic induced resistance in plants. A battery of genes is present that are involved in carbohydrate transport and metabolism, habilitating the bacterium for the assimilation of a diversity of carbon sources in the rhizosphere. Moreover, two genes, one with five copies, related to a transport system could be involved in the accumulation of glycine betaine with function of osmoprotection. Further analysis of the genome and specific genes of the 2A-2B strain will offer a better understanding of its biology, and the development of novel strategies in the biotechnological use of this bacterium.