Exploration of genetic diversity of Bacillus spp. from industrial shrimp ponds in Vietnam by multi-locus sequence typing

Abstract

Bacillus is a diverse genus consisting of more than 200 species with extensive genetic diversity. Their beneficial effects in industrial shrimp farming have been well documented. However, little is known about the biodiversity of the Bacillus spp. in this aquaculture system. Taxonomic analysis by 16S rRNA sequencing does not always allow species-level identification of Bacillus spp. In this study, 26 Bacillus isolates from two industrial Litopenaeus vannamei shrimp ponds in Bac Lieu Province, Vietnam, were analyzed for their genetic diversity by multi-locus sequence typing (MLST). A total of 22 sequence types were identified and segregated into four distinct clusters, corresponding to B. subtilis, B. velezensis, B. siamensis, and B. licheniformis. Bacillus subtilis and B. velezensis accounted for more than 73% of the Bacillus isolates. Notably, the MLST scheme exhibited high discriminatory power and might be further simplified to be a convenient method to identify species of the genus Bacillus.

Background

According to the Food and Agriculture Organization of the United Nations (FAO), aquaculture is the fastest-growing sector of food production in the world today (FAO 2018). In Vietnam, the shrimp farming area is approximately 600,000 ha, producing 300000 tons of black tiger and whiteleg shrimps per year (VASEP 2018). Although the procedure for industrial shrimp farming has been established, sustainable development of this model could be severely compromised by an increased risk of infectious diseases such as white spot syndrome virus, early mortality syndrome (EMS), and white feces syndrome.

As a result, probiotics have been increasingly employed in the form of feed supplements for shrimp farming. In Vietnam, probiotics were used in 91% of the surveyed shrimp farms (Rico et al. 2013). By definition, probiotics are live microorganisms that, when administered in adequate amounts, confer a health benefit to the host (Mack 2005). Indeed, their beneficial effects in shrimp farming have been shown in numerous studies. For instance, probiotics improve water quality, produce inhibitory compounds against pathogens, or enhance the host’s growth and immune system (Gatesoupe 1999; Gomes et al. 2009; Irianto and Austin 2002; Verschuere et al. 2000).

Bacteria belonging to Bacillus genus are often included in probiotics used in aquaculture as they are believed to confer multiple benefits to both the environment and the cultured animals (van Hai and Fotedar 2010; Zokaeifar et al. 2012). These bacteria are non-pathogenic, spore-forming, and capable of secreting compounds with antimicrobial properties (Zokaeifar et al. 2012). They have been used to promote growth and control diseases in shrimp aquaculture (Dalmin et al. 2001; Wang et al. 2005; Zokaeifar et al. 2014). However, there is a lack of knowledge on the genetic diversity of Bacillus bacteria in industrial shrimp aquaculture, which is the overall trend of shrimp farming in Vietnam.

Conventionally, culture methods or molecular techniques such as polymerase chain reaction-denaturing gradient gel electrophoresis (PCR-DGGE) (Piterina and Pembroke 2013) or 16S rRNA sequencing (Qin et al. 2016) have been used to explore the bacterial contents of aquaculture systems. However, these are time-consuming and often fail to reflect the diversity of closely related bacterial groups, particularly species of the Bacillus genus. Recently, multi-locus sequence typing (MLST), which characterizes bacterial strains using internal fragments of multiple housekeeping genes, has gained broad acceptance among epidemiologists (over 50 MLST schemes have been published and made available on the Internet at https://pubmlst.org/databases/) (Larsen et al. 2012). MLST is a standardized approach, highly unambiguous, and reproducible. Furthermore, MLST has been successfully used to study the phylogenetic diversity of the Bacillus cereus group (Sorokin et al. 2006).

In this study, we aimed to explore the genetic diversity of the Bacillus group in two industrial shrimp ponds (with and without EMS) that are frequently supplemented with probiotic products. An MLST scheme using seven housekeeping genes (glpF, ilvD, ptA, purH, pycA, rpoD, and tpiA) was applied to identify Bacillus isolates from these shrimp ponds.

Methods

Bacterial isolates

Bacillus bacteria were isolated from sediment, water, and shrimp intestine samples of two industrial whiteleg shrimp (Litopenaeus vannamei) ponds in Bac Lieu Province, Vietnam, following the procedure described by Cao et al. (2011) with some modifications. Briefly, 1 g of sample was homogenized in 100 mL of nutrient broth (NB) by Stomacher® 400 Circulator (Seward) and incubated at 80 °C for 10 min to inactivate vegetative bacteria and fungi in order to isolate Bacillus spores that withstood this heat pretreatment. The supernatant was then subjected to tenfold serial dilution before being spread onto nutrient agar (NA). After incubation at 37 °C for 24 h, individual colonies were streaked onto NA to obtain pure isolates. Upon isolation, bacterial isolates were subjected to catalase test and gram staining and positive isolates were stored in 50% glycerol at − 80 °C. A total of 26 isolates was obtained, among which 11 (sediment, n = 2; water, n = 4; intestine n = 5) were isolated from the pond that was free of EMS, while 15 (sediment, n = 8; water, n = 4; intestine n = 3) were isolated from the pond that had been affected by EMS during the last three consecutive years. Details on the origin and morphology of the isolates are presented in Table 1.

Table 1 Origin and morphology of 26 bacterial isolates used in this study

DNA extraction

DNA extraction and subsequent experiments were performed at Laboratory of Genetic Engineering, School of Biotechnology and Food Technology, Hanoi University of Science and Technology, Hanoi, Vietnam.

Total DNA of bacterial isolates was extracted following Burrell et al. (1998) with some modifications. Briefly, 2 mL of overnight LB culture was centrifuged at 10,000×g for 5 min and the supernatant was discarded. Cell pellet was then resuspended in 600 μL of Tris-EDTA (50 mM Tris pH 8.0, 5 mM EDTA). Subsequently, 50 μL of freshly prepared lysozyme (10 mg/mL) was added to the mixture and incubated at 37 °C for 2 h. A volume of 35 μL of sodium dodecyl sulfate (10% (w/v)) and 15 μL of proteinase K (10 mg/mL) was then added to the mixture, followed by another incubation step at 37 °C for 1 h. After extracting with an equal volume (700 μL) of chloroform/isoamyl alcohol (24:1, v/v), the nucleic acids from 500 μL of supernatant were precipitated by adding 50 μL of sodium acetate (3 M pH 5.2) and 1.4 mL of 100% ethanol and incubating for 1 h at room temperature. Following a centrifugation at 12,000×g for 30 min, DNA pellet was washed by 1 mL of 70% ethanol, air dried, and resuspended in 200 μL of TE (10 mM Tris pH 8.0, 1 mM EDTA) containing 10 μg/mL of RNase A. After incubating at 37 °C for 1 h to remove RNA, DNA was further purified and concentrated into a 50-μL volume using Amicon Ultra 0.5 mL 100K centrifugal filters (Millipore) following the protocols provided with the filters. DNA concentration and quality were assessed based on absorbance at 260, 280, and 230 nm using NanoDrop2000 (Thermo Fisher).

16S rRNA sequencing

16S rRNA gene of bacterial isolates was amplified by PCR using universal primers 8F (5′-AGAGTTTGATCCTGGCTCAG-3′) and 1510R (5′-GGCTACCTTGTTACGA-3′) (Ding and Yokota 2002). PCR reactions were performed with an initial denaturation at 94 °C for 3 min, followed by 30 cycles of denaturation at 94 °C for 30 s, annealing at 52 °C for 30 s, and extension at 72 °C for 1.5 min. Final extension step was performed at 72 °C for 10 min. Reaction mixtures of 50 μL contained 25 μL of GoTaq® G2 Hot Start Colorless Master Mix 2X (Promega, USA), 0.4 pmol/μL of each primer, and 10 ng of DNA template. Negative and positive (B. subtilis strain WB800N) controls were included in each PCR amplification. PCR products were purified using QIAquick PCR purification kit per the manufacturer’s specifications (QIAGEN, Germany) and sent to Macrogen (Seoul, Korea) for sequencing by Sanger method. Low-quality ends of DNA sequences were trimmed by DNA Chromatogram Explorer Lite (HeracleSoftware). DNA sequences were then BLAST searched against GenBank databases (http://www.ncbi.nlm.nih.gov) and analyzed using Bioedit (Hall 1999). MEGA X (https://www.megasoftware.net/) was used to construct the 16S phylogenetic tree using the neighbor-joining method with Kimura 2-parameter substitution model (Kikuchi 2009; Kimura 1980) and 1000 bootstrapping tests.

MLST analysis

Intragenic regions of seven housekeeping genes (glpF, ilvD, ptA, purH, pycA, rpoD, and tpiA) were selected for MLST analysis (www.pubmlst.org/bsubtilis). Primers for PCR amplification of the seven genes were designed using Primer3 software (Untergasser et al. 2012), and their sequences are shown in Table 2. PCR amplifications were performed using Promega GoTaq® G2 Hot Start Colorless Master Mix 2X as mentioned above. Reactions of 50 μL contained 25 μL of GoTaq® G2 Hot Start Colorless Master Mix 2X, 0.4 pmol/μL of each primer, and 10 ng of DNA template. One single cycling program was used for amplification of the seven genes: initial denaturation at 95 °C for 3 min, 40 cycles of denaturation (95 °C, 30 s), annealing (54 °C, 30 s), extension (72 °C, 50 s), and one final elongation step at 72 °C for 5 min. Negative and positive (B. subtilis strain WB800N) controls were included in each PCR amplification. Following amplification, PCR products were purified using QIAquick PCR purification kit or QIAquick® Gel Extraction Kit (Qiagen, Germany) per the manufacturer’s specifications and sent to Macrogen (Seoul, Korea) for sequencing.

Table 2 Primer sequences for MLST analysis

Obtained DNA sequences were trimmed at both ends to obtain regions corresponding to B. subtilis sequences available on PubMLST database (www.pubmlst.org/bsubtilis), and aligned using CLUSTALW (MEGA X). The number of polymorphic sites of each gene fragment was manually counted using the alignment outputs. Different alleles were determined on the basis of one-nucleotide difference and were assigned arbitrary numbers. For each bacterial isolate, a combination of seven alleles defined its allelic profile and sequence type (ST). Coverage of the complete coding sequences was identified using BLAST search against GenBank databases. MEGA X software was used to construct phylogenetic trees using the neighbor-joining method with Kimura 2-parameter substitution model (Akita et al. 2017; Kimura 1980) and 1000 bootstrapping tests. Sequence type analysis and recombinational tests (START) software (version 1.0.5) (http://www.mlst.net) was used to calculate G + C content and d_N/d_S value. Discrimination indices (DI) were computed as previously described (Hunter and Gaston 1988).

Results

Sequencing of 16S rRNA identified 26 Bacillus isolates

Pioneer work on prokaryotic taxonomy has recommended that identification to the species level is defined as a 16S rDNA sequence similarity of ≥ 99% with that of the type strain sequence in GenBank database (Cai et al. 2003; Stackebrandt and Ebers 2006; Benga et al. 2014). In the present study, the 16S rRNA gene fragment was amplified and sequenced using the universal primer 8F and 1510R (Ding and Yokota 2002). Approximately 1400 bp (range 1380–1421 bp) of the 16S rRNA gene sequence was successfully obtained for each isolate (Additional file 1: Table S1) with Phred scores higher than 20 (Ewing and Green 1998). These sequences were blasted against the 16S rRNA sequence database at NCBI. The results (Additional file 1: Table S1) indicated that all isolates belong to the genus Bacillus with the highest similarity scores ranging from 99.8 to 100%. However, it was not able to identify these isolates at the species level. For example, isolate BRB 2.2, BDB 1.1, BDB 11.1, BDB 3.5, BNB 1.1, BNB 1.2, BNB 5.2, BRK 5.4, BDK 2.3, BNK 2.2, BNK 2.3, BNK 7.1, and BNK 8.1 could be any species of B. amyloliquefaciens, B. velezensis, B. subtilis, or B. siamensis. The difference between the highest and second highest similarity scores was less than 0.1% for all isolates except for BRB 6.3 and BDB 6.1 (Additional file 1: Table S1).

The neighbor-joining phylogenetic tree, based on 16S rRNA sequences of the isolates and type of strains retrieved from the GenBank database, contains four clades: B. licheniformis, B. subtilis/B. tequilensis, B. amyloliquefaciens/B. siamensis, and B. velezensis (Fig. 1). From this phylogenetic tree, it is evident that the isolates BRB 6.3 and BDB 6.1 are closely related to B. licheniformis, while the isolates BNB 1.2, BNB 5.2, BNB 1.1, BRK 5.4, BDB 11.1, BNK 2.2, and BRB 2.2 are closely related to B. velezensis. Nevertheless, identification of the other isolates was inconclusive. Indeed, the low bootstrap values on the remaining part of the tree indicated that 16S rRNA sequencing is not suitable for phylogenetic analysis of all isolates at the species level (Hampl et al. 2001). This may be due to the high similarity of 16S sequences from Bacillus isolates in the present study.

Fig. 1

Neighbor-joining phylogenetic tree based on 16S rRNA sequences of 26 Bacillus isolates from EMS-free and EMS-affected shrimp ponds and representative Bacillus reference strains. Names of different clades were placed on the right-hand side. GenBank accession numbers are indicated in parentheses. Isolates from the EMS-affected pond are indicated by asterisks

All of these results clearly showed that 16S rRNA gene alone was not able to identify all Bacillus isolates at the species level. Therefore, they were subjected to genotyping by an MLST scheme that utilizes internal fragments of seven housekeeping genes.

MLST analysis

From the sequencing results, allelic and sequence profiles of the seven housekeeping genes (glpF, ilvD, ptA, purH, pycA, rpoD, and tpiA) were presented in Table 3. The lengths of analyzed fragments ranged from 384 to 470 bp, covering from 11.6 (pycA) to 55.1% (tpiA) of the complete gene sequences. Multiple sequence alignment did not show any insertions or deletions; however, SNPs were frequently observed. We found 146 (38.0%), 164 (34.9%), 105 (25.4%), 137 (34.3%), 168 (42.1%), 108 (28.1%), and 89 (21.2%) polymorphic sites for glpF, ilvD, pta, purH, pycA, rpoD, and tpiA, respectively. Moreover, for each locus, we found 11 to 19 alleles, which were counted on the basis of one-base difference. Average (G + C) content of each gene was about 49–54%. This range is similar to the (G + C) contents of the corresponding gene sequences from the B. subtilis strain 168, which is the first reference genomic data for the Bacillus genus. Average dN/dS values were much less than 1 (maximum at 0.080), indicating that the seven gene fragments are under negative selection pressure and mutations were mainly synonymous (Kryazhimskiy and Plotkin 2008). Synonymous substitutions were at least 12.5 times (1/0.080) more frequent than amino acid changes at any locus. This could be explained by the crucial functions of these housekeeping genes in Bacillus bacteria.

Table 3 Allelic profiles of the seven housekeeping genes used in MLST analysis

The discrimination indices (DI) were also computed to compare the discriminatory power of the individual genes. The lowest DI value of the seven loci was 0.908, indicating a high discriminatory power and the efficiency in differentiating the isolates in our study. glpF was scored the highest at 0.972 (18 alleles, 38.0% polymorphic sites). Interestingly, the most polymorphic fragment (pycA, 42.1% polymorphic sites) did not exhibit the highest DI (0.966). These results may allow us to further simplify the MLST scheme by using the most discriminatory loci.

After concatenation of the seven fragments, a total of 22 sequence types was distinguished among the 26 isolates. A neighbor-joining phylogenetic tree based on concatenated sequences (Fig. 2) was constructed using MEGA X software. Based on the BLAST search of the concatemer sequences, representative reference sequences were selected from GenBank databases as ingroups and outgroups. Clustering of all sequences revealed four major, non-overlapping clades, supported by the bootstrap value of 100. They corresponded to the four species of the Bacillus genus: B. velezensis, B. siamensis, B. subtilis, and B. licheniformis, respectively. We observed an uneven distribution of the isolates between these groups. Bacillus velezensis and B. subtilis clades (8 and 11 isolates, respectively) accounted for more than 73% of the total samples. Regarding the Bacillus contents in EMS-free and EMS-affected ponds, no significant difference was observed with one exception for the B. licheniformis group. Indeed, two B. licheniformis isolates were exclusively present in the EMS-affected pond. The remaining isolates from EMS-free and EMS-affected ponds were quite evenly distributed between the three clades of B. subtilis, B. velezensis, and B. siamensis.

Fig. 2

Neighbor-joining phylogenetic tree based on concatenated MLST DNA sequences of 26 Bacillus isolates from EMS-free and EMS-affected shrimp ponds and representative Bacillus reference strains. Names of different clades were placed on the right-hand side. GenBank accession numbers are indicated in parentheses. Isolates from the EMS-affected pond are indicated by asterisks

Discussion

In the present study, we described the diversity and population structure of Bacillus isolates from two industrial whiteleg shrimp ponds in Bac Lieu Province, Vietnam, by 16S rRNA sequencing and multiple-locus sequence typing. Notably, one pond was affected with EMS while the other was free of EMS. Both ponds were frequently supplemented with probiotic products.

Initially, 26 Bacillus spp. were detected by 16S rRNA sequencing. Although being useful for phylogenetic studies at the genus level, the discriminatory power at species level of the 16S method remained questionable as at least four species of the Bacillus group were identified per isolate while performing BLAST searches of the sequenced 16S fragments. This may be due to the high similarity of 16S sequences between closely related species (Stackebrandt and Goebel 1994). It has also been shown that 16S rRNA sequences of some Bacillus species are almost identical (Janda and Abbott 2007). On the other hand, the MLST scheme used in the present study allowed determination of the exact species of all the 26 isolates. Overall, all the seven genes exhibited a satisfactory discriminatory power (DI ≥ 0.908). Interestingly, the locus with the most polymorphic sites did not exhibit the highest DI (Table 2). Therefore, we suggest that using the locus with the highest discriminatory power (glpF, purH, and pycA) could be enough to differentiate bacterial isolates of B. subtilis, B. velezensis, B. siamensis, and B. licheniformis. Nevertheless, a larger population is required to evaluate this hypothesis extensively.

The neighbor-joining phylogenetic tree based on the concatenated MLST fragments showed four distinct clades corresponding to the four Bacillus species and being supported by reliable bootstrap values (higher than 80). The isolates were majorly B. subtilis and B. velezensis (73%). The dominance of B. subtilis could be due to that they are commonly used in probiotics or biocontrol agents (Buruiană et al. 2014; Farzanfar 2006). For B. velezensis, several studies have pointed out that they can act as biocontrol agents (Palazzini et al. 2016) and exhibit antimicrobial activity against fish pathogenic bacteria (Yi et al. 2018), including Vibrio parahaemolyticus, which is the leading cause of EMS in cultured shrimps. Therefore, they could have been used regularly and became widespread in industrial shrimp ponds in Vietnam. However, this is not the case for the B. licheniformis species. Despite being popular in probiotic products (Elshaghabee et al. 2017), only two isolates of this species were found in EMS-affected pond. Nevertheless, we cannot exclude the possibility that bacterial isolates identified in this study could also originate from natural sediments in the ponds. In fact, Bacillus spp. are ubiquitous and found abundantly in soil (Garbeva et al. 2003).

All Bacillus species detected in this study have been previously shown to have beneficial effects in aquaculture systems. For instance, B. subtilis and B. licheniformis are commonly used in commercialized probiotic products and their benefits have been thoroughly investigated (van Hai and Fotedar 2010; Zokaeifar et al. 2012). Several studies have also pointed out the effects of B. velezensis and B. siamensis as probiotics or biocontrol agents in industrial aquaculture farming (Buruiană et al. 2014; Meidong et al. 2017; Palazzini et al. 2016). They play a pivotal role in nutrient cycling, nutrition of the cultured animals, water quality, and disease control (Moriarty 1997).

The antagonist effects of Bacillus bacteria against V. parahaemolyticus, presumably the direct cause of EMS in shrimps, have been reported (Liu et al. 2015; Tran et al. 2013; Xu et al. 2013). However, there was no significant difference in Bacillus content between the EMS-free and EMS-affected shrimp ponds except that two B. licheniformis isolates were exclusively found in the EMS-affected pond. This preliminary result needs further research with a larger sample size to be confirmed. Of note, the ability of secreting antibacterial compounds is characteristic of a few Bacillus strains only (Azevedo et al. 1993; Liu et al. 2015). Therefore, antimicrobial activity to V. parahaemolyticus needs to be tested for each Bacillus isolates in order to determine whether there was a difference in antimicrobial profile between isolates from EMS-free and EMS-affected shrimp ponds.

Conclusions

In conclusion, we have shown that MLST is a more efficient phylogenetic tool than the 16S rRNA sequencing for identifying Bacillus species isolated from shrimp aquaculture. Using this approach, we have identified four major Bacillus species including B. subtilis, B. velezensis, B. siamensis, and B. licheniformis from EMS-free and EMS-affected industrial shrimp ponds, among which B. subtilis and B. velezensis accounted for more than 73% of the isolates. Further research will be dedicated to evaluate the antagonistic activity of the isolates against V. parahaemolyticus strains causing EMS.