Current Classification of the Bacillus pumilus Group Species, the Rubber-Pathogenic Bacteria Causing Trunk Bulges Disease in Malaysia as Assessed by MLSA and Multi rep-PCR Approaches
Article information
Abstract
Bacillus pumilus is the causal agent of trunk bulges disease affecting rubber and rubberwood quality and yield production. In this study, B. pumilus and other closely related species were included in B. pumilus group, as they shared over 99.5% similarity from 16S rRNA analysis. Multilocus sequence analysis (MLSA) of five housekeeping genes and repetitive elements-based polymerase chain reaction (rep-PCR) using REP, ERIC, and BOX primers conducted to analyze the diversity and systematic relationships of 20 isolates of B. pumilus group from four rubber tree plantations in Peninsular Malaysia (Serdang, Tanah Merah, Baling, and Rawang). Multi rep-PCR results revealed the genetic profiling among the B. pumilus group isolates, while MLSA results showed 98-100% similarity across the 20 isolates of B. pumilus group species. These 20 isolates, formerly established as B. pumilus, were found not to be grouped with B. pumilus. However, being distributed within distinctive groups of the B. pumilus group comprising of two clusters, A and B. Cluster A contained of 17 isolates close to B. altitudinis, whereas Cluster B consisted of three isolates attributed to B. safensis. This is the first MLSA and rep-PCR study on B. pumilus group, which provides an in-depth understanding of the diversity of these rubber-pathogenic isolates in Malaysia.
Malaysia is the world’s fifth-largest rubber producer that contributes more than 90% of global rubber production (Fox and Castella, 2013; Malaysian Rubber Board, 2020). The RRIM 3001 superclone rubber tree (Hevea brasiliensis) is the new high producing rubber clone introduced by Malaysian Rubber Board (MRB) to fulfill the demand for planting materials and produce a remarkable return of natural rubber production (Mokhatar et al., 2011). This clone generated a high yield of rubber production for three years of tapping and documented the topmost girth increment after six years of planting with an average of 10.6 cm/y (Nurmi-Rohayu et al., 2015). However, the RRIM 3001 superclone rubber production is under threat with a disease known as trunk bulges. The first report of the trunk bulges disease outbreak was reported in a rubber tree plantation by Mazlan et al. (2019). Interestingly, the trunk bulges problem only affects RRIM 3001 superclone rubber tree. This disease gradually reduces the quality of latex and rubberwood production. The symptoms started with the abundance of bulges of different sizes and resembled tumor-bacteriosis on the whole trunk rubber tree. Additionally, the bulges changed into cankers and produced bleeding lesions when it became worst. Bacillus pumilus was identified as the causal agent of this disease based on pathogenicity and biochemical tests and molecular characterization methods (Mazlan et al., 2019).
B. pumilus belongs to the B. substilis group (Berkeley et al., 2002; Branquinho et al., 2014c), where it contributed to a broad range of pharmaceutical and biotechnology application, including phytosanitary-based products (Espariz et al., 2016; Ficarra et al., 2016; Handtke et al., 2014; Pérez-García et al., 2011; Shah Mahmud et al., 2015) and as human and animal probiotics bacteria (Hong et al., 2005). It is widely identified as a pathogen of various plants and human diseases (Yuan and Gao, 2015). Earlier reports revealed that B. pumilus is a pathogenic bacteria causing diseases to various type of plants, including rubber (Mazlan et al., 2019), muskmelon (Song et al., 2018), ginger (Peng et al., 2013), scot pine (Kovaleva et al., 2015), ficus lacor (Hakim et al., 2015) and mango (Galal et al., 2006). Studies discovered that this bacterium shares over 99.5% of 16S rRNA similarity with other closely related species including B. altitudinis, B. stratosphericus, B. safensis, B. xiamenensis, B. aerophilus, and B. invictae (Branquinho et al., 2014a; Fritze, 2004; Lai et al., 2014; Liu et al., 2013). Later, B. invictae was reclassified as a heterotypic synonym of B. altitudinis (Liu et al., 2015a). B. pumilus group species are Gram-positive bacteria, motile, rod-shaped with ellipsoidal endospores (Mazlan et al., 2019; Satomi et al., 2006; Vettath et al., 2017). These bacteria are difficult to distinguish by phenotypic, biochemical characteristics and 16S rRNA gene sequence (Branquinho et al., 2014a; Fritze, 2004). Despite difficulty distinguishing all these species and frequently misnamed by each other, according to the PubMed data, all bacteria from marine communities are customarily placed in the B. pumilus group (Fritze, 2004). Several molecular characterization methods have been performed for bacterial differentiation, taxonomic resolution, and genetic diversity of B. pumilus group, such as single housekeeping gene sequence analysis of gyrB and rpoB regions (Branquinho et al., 2014a), multilocus sequence analysis (MLSA) using several housekeeping genes (Liu et al., 2013) and matrix-assisted laser desorption/ionization time-of-flight (Branquinho et al., 2014b).
MLSA is considered a powerful applied tool for the systematics in molecular microbes (Gevers et al., 2005) due to the drawback of 16S rRNA gene sequence as a phylogenetic marker, which is a short of resolution at the species level (Fritze, 2004; Lima-Bittencourt et al., 2007; Pontes et al., 2007). Besides, 16S rRNA application is insufficient to differentiate certain nearly related strains and species that contain high conservation, horizontal gene transfer (HGT), and multiple copies of heterogeneity of bacteria (Kitahara and Miyazaki, 2013; Liu et al., 2015b, 2017; Tian et al., 2015). Hence, to overcome these limitations, the MLSA method using protein-coding genes was suggested in parallel with 16S rRNA, since concatenating the sequence of several protein-encoding gene fragments can increase the precision and reliability of a phylogenetic scheme and provide a more robust tree topology (López-Hermoso et al., 2017; Meintanis et al., 2008; Pascual et al., 2010; Sawabe et al., 2013).
There are several factors to choose the housekeeping genes as the phylogenetic marker; uniqueness in genomes, wide distribution among bacteria, sequence divergence among related species, and phylogenetically informative size (Zeigler, 2003). Single gene sequence such as gyrB only indicates the evolution, and it may not show the “true” phylogenetic relationship. Furthermore, the concatenated aligned sequence could reduce the HGT weight (Macheras et al., 2011) and reflect the “true” relationship of bacterial taxa and serve a precise taxonomic identification between closely related strains (Glaeser and Kämpfer, 2015) and recombination of housekeeping genes (Timilsina et al., 2015). Previous studies documented successful MLSA to distinguish the closely related B. pumilus group species in the marine environment using seven housekeeping genes (gyrB, rpoB, pyrA, pyrE, aroE, mutL, and trpB) (Liu et al., 2013). Also, the MLSA has widely used for investigation of the taxonomic relationship and phylogenetic analysis of plant disease (Abidin et al., 2020; Ansari et al., 2019; Ntambo et al., 2019; Osdaghi et al., 2018; Otto et al., 2018; Oueslati et al., 2019; Suárez-Moreno et al., 2019; Waleron et al., 2019; Yahiaoui et al., 2017; Zarei et al., 2019; Zhang et al., 2018).
Repetitive sequence-based polymerase chain reaction (rep-PCR) fingerprinting, also known as PCR-based technique, is a useful tool to determine a wide range of bacteria and to compare the bacterial genome diversity of species, strains, serotypes, among others (Nurhayati et al., 2017; Rampadarath et al., 2015). It is a typing method that uses specific oligonucleotide primers corresponding to naturally occurring DNA sequences based on highly conserved repetitive DNA sequences dispersed throughout the genome of diverse bacterial species (da Silva and Valicente, 2013). Numerous amplicons of distinctive electrophoresis patterns establishing the DNA fingerprint-specific pattern for an individual bacterial strain (Mohapatra et al., 2007; Rademaker and de Brujin, 1997).
Most rep-PCR DNA fingerprinting studies are based on three repetitive elements, which are repetitive extragenic palindromic (REP) elements (REP-PCR) with 35-40 bp sequences, enterobacterial repetitive intergenic consensus (ERIC) elements (ERIC-PCR) and box A, B, and C subunits (BOX) element (BOX-PCR) with 124-127 bp sequences and 154 bp sequences, respectively (Louws et al., 1994, 1999; Pasanen et al., 2014; Rademaker et al., 2004). The rep-PCR method has frequently distinguished the closely related species within genus Bacillus such as B. cereus, B. thuriengiensis, and B. anthracis (Cherif et al., 2003, 2007; da Silva and Valicente, 2013). It is a simple, reliable, reproducible and highly sensitive method for analyzing the distribution of repetitive DNA sequence in several prokaryotic genomes (Amoupour et al., 2019; Louws et al., 1999; Masanto et al., 2019; Rademaker and de Brujin, 1997; Versalovic et al., 1991).
A collection of pathogenic isolates, previously known as B. pumilus, infecting RRIM 3001 superclone rubber tree was analyzed in this present work (1) to determine the genetic relationship of B. pumilus group isolates associated with trunk bulges disease of RRIM 3001 superclone rubber tree from different geographical areas in Peninsular Malaysia by MLSA and (2) to elucidate the genetic diversity of B. pumilus isolates associated with trunk bulges disease using molecular profiling via the rep-PCR fingerprinting method. These results represent the first MLSA and multi rep-PCR studies on B. pumilus group species of rubber-pathogenic isolates that may pave the way in the current taxonomic classification of these pathogens. They are crucial for developing B. pumilus group species’ management strategies in rubber plantations since rubber tree has been recognized as an essential commodity crop with a high economic value in Malaysia.
Materials and Methods
Isolates collection and identification
A total of 20 B. pumilus isolates used for phylogeny study and rep-PCR method; eight isolates from Serdang in Selangor, five isolates from Baling in Kedah, five isolates from Rawang in Selangor, and two isolates from Tanah Merah in Kelantan. These isolates were collected from the different outbreak trunk bulges rubber tree plantations in Peninsular Malaysia and earlier confirmed as B. pumilus based on phenotypic and molecular characterization (Mazlan et al., 2019). The bacterial isolates were stored in glycerol at ‒20°C. Detailed information of 20 isolates of B. pumilus group is listed in Table 1.
DNA extraction
Bacillus pumilus group isolates were grown in nutrient broth for 24 to 48 h at 28°C. Twenty bacterial isolates of B. pumilus group extracted using a commercial genomic DNA isolation kit (Presto Mini gDNA Bacteria Kit, Geneaid Biotech Ltd., New Taipei City, Taiwan) following the protocol provided.
PCR amplification and sequencing of 16S rRNA and housekeeping genes
PCR was performed in 25 µl reaction mixture, containing 3 µl of genomic DNA template, 12.5 μl of 2× DreamTaq Red PCR MasterMix (Thermo Scientific Inc., Waltham, MA, USA), 0.5 µl of each primer, and 8.5 µl of sterile distilled water. PCR amplification was performed using an ‘iCycler’ Thermal Cycler (BioRad Laboratories Inc., Hercules, CA, USA) to amplify the 16S rRNA gene using universal pair primers 27F (5′-AGA GTT TGA TCC TGG CTC AG-3′) and 1492R (5′-ACGGCTACCTTGTTACGACT-3′), and five housekeeping genes primer sets (gyrB, pyrE, aroE, rpoB, and trpB) of B. pumilus group (Liu et al., 2013) (Table 2).
List of primers used in the MLSA study
Each PCR product was separated by electrophoresis on 1% agarose gel at 80 V (Power supply model 1000/500, Bio-Rad Laboratories Inc.) for 35 min using Mini Sub DNA Cell (Bio-Rad Laboratories Inc.) and imaged by Alpha Imager System (Alpha-Innotech, Siber-Hegner, UK). All amplified PCR products were sent for sequencing (MyTACG Bioscience Enterprise, Selangor, Malaysia). Sequences were deposited in the NCBI GenBank database for accession numbers (Table 3).
GenBank accession numbers of all 20 isolates of Bacillus pumilus group used in this study
Phylogenetic analysis of 16S rRNA, single and concatenated housekeeping genes
The determined sequences of the 16S rRNA gene and five housekeeping genes were analyzed against sequences in the NCBI database using BLASTn (Altschul et al., 1990). Both were aligned using the ClustalW algorithm with manual adjustments implemented in MEGA 7.0 software (Kumar et al., 2016). The genetic distance and sequence similarity were calculated using Kimura’s 2-parameter model (Kimura, 1980) with MEGA 7.0 software. The polymorphic sites (S) were calculated using DnaSP 5.0 software (Librado and Rozas, 2009). Phylogenetic trees were constructed using Maximum-likelihood analysis and the model of evolution was determined with the MODELTEST tab in MEGA 7.0 software (Hall, 2013). The strength of the resulting trees‘ internal branches was evaluated by bootstraps analysis with 1000 bootstrap replication. Four strains types, B. altitudinis MCCC1A01287T, B. safensis MCCC1A00456T, B. pumilus MCCC1A00439T and MCCC1A06996T were included in the phylogeny study. However, B. stratosphericus, B. aerophilus, B. xiamenensis, and B. invictae were not available in the public collection (Liu et al., 2013). B. cereus ATCC 14579T (GenBank accession no. AE016877) was used as an outgroup.
REP-PCR, ERIC-PCR, BOX-PCR and rep-PCR amplification
DNA of the 20 B. pumilus group isolates was subjected to rep-PCR using three primer sets (De Bruijn, 1992; Louws et al., 1994; Thwaites et al., 1999): REP (REP1R-I, 5′-IIIICGYCGICATCMGGC-3′; REP2-I, 5′-ICGICTTATCIGCCGGTAC-3′; ERIC (ERIC1R, 5′-ATGTAAGCTCCTGGGATTCAC-3′; ERIC2, 5′-AAGTAAGTGACTGGGGTGAGC-3′) and BOX (BOXA1R, 5′-CTACGGCAAGGCGACGCTGACG-3′). The reaction mixture (25 μl) consisted of 5 μl of template DNA, 12.5 μl of Taq DNA polymerase, 5.5 μl of sterilized distilled water, and 1 μl of each primer. PCR reactions were carried out with cycling conditions as follows: 5 min denaturation at 94°C followed by 30 cycles of 94°C for 1 min, annealing at 43.5°C, 47.5°C, or 52.4°C for 1 min for REP, ERIC, and BOX primer sets, respectively, and 8 min extension at 65°C, with a final extension of 15 min at 65°C. PCR fragments were electrophoretically separated on 2% agarose gel with FloroSafe DNA stain at 75 V (Power supply model 1000/500, Bio-Rad Laboratories Inc.) for 45 min using Mini Sub DNA Cell (Bio-Rad Laboratories Inc.), then imaged by Alpha Imager System (Alpha-Innotech).
Each primer set was repeated at least three times to confirm amplification consistency, and only clear and reproducible banding patterns were used for phylogenetic analysis. The clear, unambiguous, and reproducible monomorphic and polymorphic bands were scored and evaluated by recapitulating into 0-1 table (0, no band; 1, band present) and putting the band arrangement in columns and isolate number in rows using Excel (Microsoft, Redmond, WA, USA). Diversity analyses were generated in numerical taxonomy and multivariate analysis system (NTSYSpc version 2.21q) (numerical taxonomy system, Applied Biostatistics, Port Jefferson, NY, USA). Similarity among 20 B. pumilus group species was established as matrices of genetic similarity compiled using the SIMQUAL function Jaccard‘s coefficients. Dendrograms representing the genetic relationship among all B. pumilus group species were generated from the similarity matrices by applying the unweighted pair-group arithmetic mean method (UPGMA) (cluster analysis) with the SAHN function system.
Results
Phylogenetic analysis of 16S rRNA and five single housekeeping genes
PCR amplification of all 20 isolates using 16S rRNA gene analysis, each produced ~900 bp amplicon. 16S rRNA gene analysis was unable to distinguish B. safensis, B. altitudinis and B. pumilus isolates from each other. All 20 isolates, formerly identified as B. pumilus, being distributed within distinctive groups of the B. pumilus group (Fig. 1). Seventeen isolates were close to the type strain B. altitudinis with a bootstrap support value of 91% and three isolates were close to the type strains B. safensis and B. pumilus, with a bootstrap support value of 98%.
Phylogenetic tree based on the 16S rRNA gene sequence comparison, showing the relationship of Bacillus pumilus group species; B. altitudinis, B. safensis, and B. pumilus using a maximum-likelihood method with MEGA 7.0. Bacillus cereus ATCC 14579 was used as an outgroup. All isolates used in this study was signified using a purple colored rectangle. Scale bar = 0.005 substitutions per nucleotide position.
Phylogenetic affiliation based on five housekeeping genes (gyrB, pyrE, aroE, rpoB, and trpB) successfully distinguished the isolates and obtained a precise overview of the phylogenetic position of the bacterial isolates. All 20 B. pumilus group isolates delineated into two distinct phylogenetic clusters, B. altitudinis and B. safensis. None of the collected 20 isolates was grouped into B. pumilus MCCC1A00439T and MCCC1A06996T strains. Specifically, these 20 isolates of B. pumilus group were divided into two clusters, A and B. Cluster A contained 17 isolates (SD1, SD2, SD11.1, SD12.4, SD12.1, R2.3, R2.5, R7.2, R7.1, R3.1, KD14.4, KD15.5, KD15.2, KD3.1, KD14.1, KEL1.2, and KEL1.1) close to type strain B. altitudinis, while Cluster B consisted of 3 isolates (SD11.3, SD12.6, and SD12.9) attributed to the type strain B. safensis. All individual phylogenetic trees showed almost similar topology structure and revealed 89-100% similarities among the 20 B. pumilus group bacterial isolates.
As shown in previous studies, gyrB gene analyses exhibit a high-resolution power because this gene produced a faster evolution rate inferred from 16S rRNA gene sequences. Similarly, in this study, the phylogenetic tree of the gyrB gene (Fig. 2) divided 20 isolates of B. pumilus group into two clusters, Cluster A and B. Cluster A composed of 17 bacterial isolates (SD1, SD2, SD11.1, SD12.4, SD12.1, R2.3, R2.5, R7.2, R7.1, R3.1, KD14.4, KD15.5, KD15.2, KD3.1, KD14.1, KEL1.2, and KEL1.1) attributed to B. altitudinis with 91% bootstrap value and 3 isolates (SD11.3, SD12.6, and SD12.9) clustered together with B. safensis with 98% bootstrap value in Cluster B.
Phylogenetic tree based on the gyrB gene sequence comparison of all 20 isolates of Bacillus pumilus group using a maximum-likelihood method with MEGA 7.0. All isolates was classified into Cluster A and Cluster B. Bacillus cereus ATCC 14579 was used as an outgroup. Scale bar = 0.10 substitutions per nucleotide position.
There were slight differences observed in some topologies of the five phylogenetic trees. In gyrB and trpB phylogenetic trees, B. pumilus is closer to B. altitudinis, in rpoB and pyrE phylogenetic trees, B. altitudinis closer to B. safensis. Unlike other trees, aroE phylogenetic tree showed that B. altitudinis closer to B. safensis. Other minor differences were also observed in the trees, as shown in the Supplementary Figs. 1-5.
Characteristics of 16S rRNA and five housekeeping genes
All 20 isolates of B. pumilus group were analyzed to discriminate the closely related species. Characteristics of each housekeeping gene, including the gene length, polymorphic sites, the mean G + C content and the genetic distance are listed in Table 4.
Characteristics of 16S rRNA gene and five housekeeping genes from all 20 isolates of Bacillus pumilus group used in this study
The genetic distance of all housekeeping genes exhibited high-resolution power compared to the 16S rRNA gene. Among all the housekeeping genes used, aroE showed the highest resolution power with 17.69% polymorphic sites and the most extensive genetic distance range of 0.000-0.181 (mean, 0.059), while rpoB showed the lowest resolution power with 4.617% polymorphic site and a genetic distance range of 0.000-0.035 (mean, 0.020).
Phylogenetic analysis based on the concatenated housekeeping genes
All five housekeeping genes were concatenated in the order of gyrB-pyrE-aroE-rpoB-trpB to reexamine the phylogeny of the 20 B. pumilus group isolates. This new phylogenetic tree showed a similar topology as the tree described above based on the five single genes but was more elaborated and stable (Fig. 3). Specifically, Cluster A consisted of 17 isolates (SD1, SD2, SD11.1, SD12.4, SD12.1, R2.3, R2.5, R7.2, R7.1, R3.1, KD14.4, KD15.5, KD15.2, KD3.1, KD14.1, KEL1.2, and KEL1.1) belonging to B. altitudinis with 98% bootstrap value. Cluster B consisted of three isolates (SD11.3, SD12.6, and SD12.9) belonging to B. safensis with 100% bootstrap value.
Phylogenetic tree based on five housekeeping genes (gyrB-pyrE-aroE-rpoB-trpB) sequence comparison of all 20 isolates of Bacillus pumilus group using a maximum-likelihood method. All isolates was classified into Cluster A consisted of 17 isolates belonging to B. altitudinis and Cluster B consisted of three isolates belonging to B. safensis. Bacillus cereus ATCC 14579 was used as an outgroup. Scale bar = 0.10 substitutions per nucleotide position.
REP, ERIC, BOX-PCR, and rep-PCR amplification
All three primer sets generated specific and clear fingerprinting patterns to amplify genomic DNA for B. altitudinis and B. safensis isolates (Fig. 4). REP, ERIC, and BOX primers generated 202, 200, and 49 bands, respectively. In rep-PCR analyses, a 1,400 bp amplified fragment was relatively specific to B. altitudinis isolates, while B. safensis isolates had a relatively uniform rep-PCR fingerprinting pattern with significant bands at about 100 bp, 500 bp, and 1,500 bp sizes. ERIC-PCR analyses revealed that 250 bp and 500 bp fragments were amplified, which relatively specific to B. altitudinis, while B. safensis isolates had a relatively uniform rep-PCR fingerprinting pattern with major bands of about 350 bp and 550 bp sizes. BOX-PCR analyses later disclosed no fingerprinting pattern on the B. altitudinis isolates but produced a 1000 bp amplicon specific for all B. safensis isolates.
Multi rep-PCR fingerprinting analyses of all 20 isolates of Bacillus pumilus group isolates using REP-PCR (A), ERIC-PCR (B), and BOX-PCR (C). rep-PCR, repetitive elements-based polymerase chain reaction; REP, repetitive extragenic palindromic; ERIC, enterobacterial repetitive intergenic consensus; BOX, box A, B, and C subunits.
Dendrogram of combined data of REP-, ERIC-, and BOX-PCR amplification of multi rep-PCR data revealed that all 20 isolates were grouped into two groups. Cluster A consisted of 17 total isolates from Baling, Tanah Merah, Rawang and five isolates from Serdang (SD1, SD2, SD11.1, SD12.4, and SD12.6). Cluster B comprised only three isolates (SD11.3, SD12.1, and SD12.9) from Serdang. The percentage of similarity coefficient for Cluster A was ~34%, and ~90% for Cluster B (Fig. 5).
Dendrogram showing diversity of all 20 isolates of Bacillus pumilus group based on multi rep-PCR (REP-PCR, ERIC-PCR, and BOX-PCR) analyses using the UPGMA clustering method. All isolates was classified into Cluster A and Cluster B. rep-PCR, repetitive elements-based polymerase chain reaction; REP, repetitive extragenic palindromic; ERIC, enterobacterial repetitive intergenic consensus; BOX, box A, B, and C subunits; UPGMA, unweighted pair-group arithmetic mean method.
Discussion
This research represents the first report on the genetic diversity of B. pumilus group isolates associated with trunk bulges disease of RRIM 3001 superclone rubber tree (H. brasiliensis) in Malaysia, using both MLSA and rep-PCR fingerprinting methods. B. pumilus group consisted of seven species that share over 99.5% of 16S rRNA gene identity and are frequently misnamed. Hence, developing a reliable, rapid and affordable alternative technique is crucial for accurate taxonomic B. pumilus group species (Branquinho et al., 2014a).
Based on these findings, all individual phylogenetic trees (gyrB, pyrE, aroE, rpoB, and trpB) and concatenated sequence genes tree confirmed the allocation of all 20 isolates to the B. pumilus group species. All single and concatenated phylogenetic trees exhibited almost identical topology, although some variations were found among the different genes analyzed. Phylogenetic tree analysis from each gene revealed that these 20 isolates that belonged to the B. pumilus group were clustered into two groups with around 89-100% similarity and differentiated into 3 group species; B. pumilus, B. altitudinis, and B. safensis. These results corroborated with previous findings on proteinencoding genes‘ capability to clarify the discrimination between closely related species and generate robust tree topology compared to ribosomal gene analysis, which shows low-resolution power (Azevedo et al., 2015; Brady et al., 2008; Sawabe et al., 2013). Ng (2020) documented that ribosomal proteins that hold partial phylogenetic significance constitute the second group. The ribosomal proteins could reproduce significant branches of the 16S rRNA phylogenetic tree but had difficulty differentiating B. licheniformis and B. pumilus at the sequence level.
Among the five housekeeping genes, aroE possessed a relatively high differentiation power based on the most extensive genetic distance range at 0.000-0.181 (mean, 0.059) and the highest polymorphic sites of 17.69%, followed by pyrE gene with 16.84% polymorphic site and 0.000-0.180 genetic distance range (mean, 0.053). Both aroE and pyrE genes possessed the highest resolution power compared to gyrB, rpoB and trpB genes. However, gyrB gene has widely been applied for taxonomic resolution of closely related species and is commonly recommended for Bacillus species classification due to the invigorating popularity of gyrB gene in the GenBank database (Adékambi et al., 2008; Fritze, 2004; Konstantinidis et al., 2006). Based on the promising results, we suggest that the aroE, pyrE and gyrB can potentially be utilized as standard markers to distinguish the closely related species of B. pumilus group.
Concatenations of multiple loci in phylogenetic tree defeat many single-gene analyse’ limitations and define solid and robust relationships for bacterial species classification (Hanage et al., 2005; Papke et al., 2011). In this study, the concatenated phylogenetic revealed that all isolates grouped into B. altitudinis and B. safensis species clusters with 98% and 100% bootstrap values (Fig. 3). Cluster A consisted of 17 isolates including SD1, SD2, SD11.1, SD12.4, SD12.1, R2.3, R2.5, R7.2, R7.1, R3.1, KD14.4, KD15.5, KD15.2, KD3.1, KD14.1, KEL1.2, and KEL1.1 were grouped with B. altitudinis reference strain, while three isolates, SD11.3, SD12.6, and SD12.9) were clustered into B. safensis reference strain in Cluster B. Unexpectedly, all 20 isolates used to be formerly identified as B. pumilus, actually belonged to the species of B. altitudinis and B. safensis. These isolates collected from four collection areas not clustered together, showing high variation among isolates from similar geographical locations. However, isolates from different geographical locations exhibited less genetic variations. The results provide evidence that phylogenetic analyses using MEGA 7.0 successfully formed a robust and distinct cluster between all 20 isolates with B. pumilus MCCC1A00439T and MCCC1A06996T reference strains, B. safensis MCCC1A04526T reference strain and B. altitudinis MCCC1A01287T reference strain. Similar results obtained from previous MLSA studies conducted on distribution of B. pumilus group members in marine environments using phylogenetic trees based on multiple housekeeping genes (gyrB, rpoB, aroE, pyrE, mutL, pyrA, and trpB), which successfully distinguish between closely related species of B. pumilus group (Liu et al., 2013).
Many reports disclosed that some B. pumilus group species are pathogenic bacteria, especially B. pumilus and B. altitudinis. B. pumilus is a pathogenic bacteria-causing disease to various crops, including fruit rot of muskmelon (Song et al., 2018), rhizome rot of ginger (Peng et al., 2013), soft rot of scot pine seedling (Kovaleva et al., 2015), fruit rot of ficus lacor (Hakim et al., 2015) and leaf blight of mango tree (Galal et al., 2006). Bacillus altitudinis was reported as the new causal agent of soft rot pathogen for both apple and pear fruits, in which most of its virulent strains causing high soft rot severity on apple and pear cultivars based on the disease severity index, phenotypic test and molecular characterization (Elbanna et al., 2014). On the contrary, B. safensis is continuously used as a plant growth-promoting rhizobacteria (Khan et al., 2017; Wu et al., 2019). Multiple genome comparisons signified that some B. safensis isolates have mistakenly been identified as B. pumilus, especially when extensive molecular analyses were not considered (Tirumalai et al., 2018). Agbobatinkpo et al. (2013) reported that B. safensis shared 90.2% gyrA sequence similarity with B. pumilus, which is nearly the same as a result (91.2% gyrB sequence similarity) obtained by Satomi et al. (2006). A recent phylogenetic study on B. pumilus and B. safensis strains FO-36b and MERTA revealed that both strains were clustered together in a distinct group of B. safensis (Espariz et al., 2016). Branquinho et al. (2014b) suggested that ribosomal and spore proteins constituted most B. pumilus and B. safensis biomarkers, whose fingerprinting by matrix assisted laser desorption ionization-time of flight mass spectrometry and other MSbased techniques can be used for rapid and accurate identification of B. safensis.
Glaeser and Kämpfer (2015) documented that MLSA analysis should be carried out on more than four housekeeping genes for a more stable, precise, and accurate topology tree. Our study discovered that MLSA was insufficient to clear up these issues and need fingerprinting tools to completely differentiate closely related species of B. pumilus group as the housekeeping genes only occupy 0.1-0.2% of the genome (Liu et al., 2013). Hence, the repPCR method using REP, ERIC, and BOX primers proved to generate more accurate information due to the higher variability of its genetic region than other genomic regions for analyzing genetic diversity and the relationship between closely related isolates obtained in this study. The rep-PCR method has been considered an excellent discriminatory tool and more reproducible than other fingerprinting methods for analyzing the genetic diversity among the isolates of Bacillus sp. (Cherif et al., 2003, 2007; da Silva and Valicente, 2013). The presence of similar banding patterns as observed from at least three PCR replications demonstrate the reproducibility of this technique and its suitability for use in B. pumilus group genetic studies.
In rep-PCR analyses, all three primer sets generated specific and clear fingerprinting patterns for B. altitudinis and B. safensis isolates (Fig. 4). Both of them exhibited their fingerprinting pattern. Interestingly, distinct fingerprinting patterns were shown among the two species of B. pumilus group isolates collected from different geographical locations in Peninsular Malaysia. Based on the primers results, REP-PCR is the most suitable method for species clustering and separation since REP primers displayed the most apparent fingerprinting pattern, pursued by ERIC primers and BOX primers.
The dendrograms generated by REP-, ERIC-, BOXPCRs showed that both species separated well. Subsequently, all data from the three primers (REP, ERIC, and BOX) were combined in a single similarity matrix for groupings’ validity, collectively known as multi rep-PCR analysis. The multi rep-PCR analysis is sensitive enough to characterize the relationships and assess intraspecific diversity. Compared to independent rep-PCR experiments, multi rep-PCR generates a higher level of discrimination, whereas isolates of the two species remain separate in the corresponding dendrogram. Based on our findings, the two species of B. pumilus group appeared well-separated. Our findings discovered the similarity of genetic profiles within clusters of the B. pumilus group isolates, thus recommending a different status of these two species since all B. altitudinis isolates were clustered into Cluster A, while all B. safensis isolates were clustered into Cluster B. For REP-, ERIC-, BOX- and multi rep-PCRs topology analysis, Cluster A represented all B. altitudinis isolates from Kedah, Kelantan, Rawang, Selangor, and five isolates from Serdang, Selangor, while Cluster B represented three B. safensis isolates from Serdang, Selangor. Three isolates of B. altitudinis (SD1, SD2, and SD12.4) from Serdang, Selangor and three isolates of B. altitudinis (KD14.4, KD15.5, and KD14.1) from Baling, Kedah exhibited similar genotype pattern for all three primer sets. Meanwhile, two isolates of B. altitudinis (R7.2 and R7.1) from Rawang, Selangor exhibited similar genotype patterns for all three primer sets. This distribution is in agreement with previous works based on MLSA results obtained from this study.
The next appealing issue raised from the multi-rep-PCR analysis is the intraspecific diversity and isolate distribution within the species. We have disclosed the similarity of genetic profiles within clusters of the B. pumilus group isolates based on the dendrogram of the combined data from multi rep-PCR. The similarity was observed among isolates in Cluster A (all B. altitudinis isolates from Kedah, Kelantan, Rawang, Selangor and five isolates of Serdang, Selangor isolates; ~35%) and Cluster B (three B. safensis isolates from Serdang, Selangor; ~95%) (Fig. 5) from different geographical locations. The results signified that the isolates conceivably were possibly lineages of a single virulent isolate and that transmission with planting stock or infected seeds was a likely means of spread across Peninsular Malaysia. All 20 bacterial isolates were clustered together, indicating that the isolates were genetically similar and shared close phylogenetic relatedness. All 20 isolates were probably commenced and derived from a single emergence a long time ago (Moretti et al., 2017). Almost all 17 B. altitudinis isolates had similar band patterns and were grouped in a separate dendrogram branch. All B. altitudinis isolates were grouped in Cluster A together, revealing that these isolates derived from the same geographical regions relatively homogeneous. The overall distribution of B. safensis type isolates and their apparent separation in Cluster B of the dendrograms validated the species’ high phylogenetic homogeneity. The rep-PCR fingerprinting has proven as an excellent tool to characterize and discriminate the B. pumilus group isolates at the genomic level. Kumar et al. (2014) supported this study and Patil et al. (2010) that these polyphasic genotypic fingerprinting techniques are excellent and reliable tools discriminating Bacillus isolates as a separate group.
The rep-PCR method provided more discriminatory power than the MLSA method, based on the fingerprinting profile, to differentiate between all 20 B. pumilus isolates. Although both of the methods, MLSA and rep-PCR, provided greater insight than one method alone, there are distinct advantages in using MLSA to determine the genetic relationships. MLSA generates a discrete data set based on the nucleotide sequences of known genes and allows for accurate calculate of genetic distances compared to repPCR. The results are portable and additional sequences can be added to the database as they become available.
These findings may add new knowledge on the distribution and the genetic diversity of B. pumilus group isolates associated with trunk bulges disease of RRIM 3001 superclone rubber tree in Peninsular Malaysia. This study clearly showed that the isolates of B. pumilus group from different geographical regions and locations in Peninsular Malaysia were distinct among isolates from a similar geographical location. However, between isolates from different geographical locations, the variations were much less. The MLSA and multi rep-PCR methods were successful in distinguishing the B. pumilus group species from each other. Generating distribution and diversity maps would provide valuable information to design disease control strategies to limit B. pumilus group species‘ introduction into new regions or host plants. Moreover, we have provided crucial information to plant breeders seeking to incorporate durable tolerance or resistance into commercial cultivar by profiling the taxonomic diversity among the B. pumilus group isolates.
Notes
Conflicts of Interest
No potential conflict of interest relevant to this article was reported.
Acknowledgements
We acknowledge Universiti Putra Malaysia (GPB-9552900) and the Fundamental Research Grant Scheme from the Ministry of Education Malaysia (FRGS-5540090) for financial supports. Ainur Ainiah Azman Husni is a recipient of the Graduate Research Fellowship from Universiti Putra Malaysia.
Electronic Supplementary Material
Supplementary materials are available at Plant Pathology Journal website (http://www.ppjonline.org/).
Phylogenetic tree based on the gyrB gene sequence comparison of all 20 isolates of Bacillus pumilus group using a maximum-likelihood method with MEGA 7.0. All isolates were classified into Cluster A and Cluster B. Bacillus cereus ATCC 14579 was used as an outgroup. Scale bar = 0.10 substitutions per nucleotide position.
Phylogenetic tree based on the pyrE gene sequence comparison of all 20 isolates of Bacillus pumilus group using a maximum-likelihood method with MEGA 7.0. All isolates were classified into Cluster A and Cluster B. Bacillus cereus ATCC 14579 was used as an outgroup. Scale bar = 0.10 substitutions per nucleotide position.
Phylogenetic tree based on the aroE gene sequence comparison of all 20 isolates of Bacillus pumilus group using a maximum-likelihood method with MEGA 7.0. All isolates was classified into Cluster A and Cluster B. Bacillus cereus ATCC 14579 was used as an outgroup. Scale bar = 0.10 substitutions per nucleotide position.
Phylogenetic tree based on the rpoB gene sequence comparison of all 20 isolates of Bacillus pumilus group using a maximum-likelihood method with MEGA 7.0. All isolates were classified into Cluster A and Cluster B. Bacillus cereus ATCC 14579 was used as an outgroup. Scale bar = 0.10 substitutions per nucleotide position.
Phylogenetic tree based on the trpB gene sequence comparison of all 20 isolates of Bacillus pumilus group using a maximum-likelihood method with MEGA 7.0. All isolates were classified into Cluster A and Cluster B. Bacillus cereus ATCC 14579 was used as an outgroup. Scale bar = 0.10 substitutions per nucleotide position.