De novo assembly and comparative genome analysis for polyhydroxyalkanoates-producing Bacillus sp. BNPI-92 strain

Abstract

Background: Certain Bacillus species play a vital role in polyhydroxyalkanoate (PHA) production. However, most of these isolates did not properly identify to species level when scientifically had been reported.

Results: From NGS analysis, 5719 genes were predicted in the de novo genome assembly. Based on genome annotation using RAST server, 5,527,513 bp sequences were predicted with 5679 bp number of protein-coding sequence. Its genome sequence contains 35.1% and 156 GC content and contigs, respectively. In RAST server analysis, subsystem (43%) and non-subsystem coverage (57%) were generated. Ortho Venn comparative genome analysis indicated that Bacillus sp. BNPI-92 shared 2930 gene cluster (core gene) with B. cereus ATCC 14579 ^T (AE016877), B. paranthracis Mn5T (MACE01000012), B. thuringiensis ATCC 10792 T (ACNF01000156), and B. antrics Amen T (AE016879) strains. For our strain, the maximum gene cluster (190) was shared with B. cereus ATCC 14579 ^T (AE016877). For Ortho Venn pair wise analysis, the maximum overlapping gene clusters thresholds have been detected between Bacillus s p.BNPI-92 and Ba. cereus ATCC 14579 ^T (5414). Average nucleotide identity (ANI) such as OriginalANI and OrthoANI, in silicon digital DND-DNA hybridization (isDDH), Type (Strain) Genome Server (TYGS), and Genome-Genome Distance Calculator (GGDC) were more essentially related Bacillus sp. BNPI-92 with B. cereus ATCC 14579 ^T strain. Therefore, based on the combination of RAST annotation, OrthoVenn server, ANI and isDDH result Bacillus sp.BNPI-92 strain was strongly confirmed to be a B. cereus type strain. It was designated as B. cereus BNPI-92 strain. In B. cereus BNPI-92 strain whole genome sequence, PHA biosynthesis encoding genes such as phaP, phaQ, phaR (PHA synthesis repressor phaR gene sequence), phaB/phbB, and phaC were predicted on the same operon. These gene clusters were designated as phaPQRBC. However, phaA was located on other operons.

Conclusions: This newly obtained isolate was found to be new a strain based on comparative genomic analysis and it was also observed as a potential candidate for PHA biosynthesis.

Keywords: ANI; Bacillus species; Gene; Genome DNA; Metabolic pathway.

Fig. 1

Graphical representation of annotated genome assembly of BNPI-92 based on RAST server. Genome assembly (genes) connected to subsystems and their distribution in different categories. a Subsystem coverage, b subsystem category of percentage distribution, and c subsystem feature count that expandable down to the specific gene with their respective role (see secondary metabolism). This online tool is available at http://rast.theseed.org/FIG/rast.cgi

Fig. 2

Graphic circular genome mapping for Bacillus sp. BNPI-92 that produced PHA polymers and obtained from area of plastic wastes accumulation. This genome map is visualized by CGview.ca that available at http://stothard.afns.ualberta.ca/cgview_server/

Fig. 3

Ortho Venn analysis using online tools that available at available via https://orthovenn2.bioinfotoolkits.net/task/create. a The occurrence table contains multiple groups of gene cluster (the pattern to the left which indicates the species are in the clusters) such as cluster count (number of gene clusters shared between species) and protein count (number of protein members in the shared cluster for these strains). Row indicates orthologous gene cluster for multiple species that summarized as a cell graph and column indicates different closely related bacterial species. The occurrence table with deep purple color bar represented the pattern of shared multiple orthologous gene cluster among Bacillus sp. BNPI-92 and other closely related bacterial strains whereas gray color bars indicate the absence of gene cluster in these strains (Fig. 5a). b OrthoVenn diagram graphic tools used for comparing a protein sequence of PHA-producing BNPI-92 strain with other five closely related strains. c Similarity matrix for pairwise protein sequence comparison for heatmap that shows the orthologs cluster between Bacillus sp. BNPI-92 and other closely related strains of protein sequence. This heatmap [43] was computed into ClustVis server that is available at https://biit.cs.ut.ee/clustvis/

Fig. 4

a Top 10 terms in biological process category from GO annotation of BNPI-92 strain. b Top 10 terms in molecular function category from GO annotation of BNPI-92 strain. c Top cellular component category from GO annotation of BNPI-92 strain

Fig. 5

a Phylogenetic tree and evolutionary relationships of taxa for Bacillus sp. BNPI-92 and the other closely related strains. The evolutionary history was inferred using the neighbor-joining method [45]. The optimal tree with the sum of branch length = 0.18016536 is shown. Evolutionary distances were calculated by neighbor-joining and based on 1000 bootstrap replication of confidence values (percentage of 1000 replication). Bar, 0.05 substitutions per nucleotide position [46]. The evolutionary distances were computed by using p-distance method [47]. The analysis involved six nucleotide sequences. All positions containing gaps and missing data were eliminated. Finally, evolutionary analyses were performed with MEGA7.0.9 software [48]. Micrococcus luteus DSM 20030 T (AJ536198.1) was designated and used as outgroup in the analyses; other related sequences were obtained from EzTaxon-e server and annotated using RAST server before tree construction b UPGMA dendrograms heatmap for OrthoANI [17], c UPGMA dendrograms heatmap for OriginalANI [17] and d UPGMA dendrograms heatmap for genome to genome distance calculator (GGDC) [17]. e In silicon DDH ( is DDH) for Bacillus sp. BNPI-92 and closely related strains. In silicon DDH analysis was performed using online tool that is available at https://tygs.dsmz.de/and Fig. 8e constructed using GraphPad Software. f Average nucleotide identity (ANI) (%) between pairs of Bacillus sp. BNPI and other five strains. ANI values of ≥ 96% and is DDH values of ≥ 70% consistently grouped genomes originating from strains of the same species together. It was performed according to Goris et al . [31] and using online tools that are available at http://enve-omics.ce.gatech.edu/ani/

Fig. 6

GBDP phylogeny based on genome data and TYGS result for PHA-producing bacterial isolate data set. Tree inferred with Fast ME 2.1.6.1 [37] from GBDP distances calculated from genome sequences of Bacillus sp. BNPI-92 and other closely related Bacillus strains. The branch lengths are scaled in terms of GBDP distance formula d₅. The numbers above branches are GBDP pseudo-bootstrap support values > 60% from 100 replications, with average branch support of 99.9%. The tree was rooted at the midpoint [38]. A labeled and colored box are annotated by affiliation and corresponding to (1) isDDH species clusters, (2) isDDH sub-species clusters, (3) percent GC (34.8–35.3%), (4) delta statistics (0.1–0.2), (5) genome size (4,614,627–6,234,842 bp), (6) protein count (5255–6243), and (7) user strain (query sequence)

Fig. 7

Annotation overview and schematic representation of secondary metabolite gene clusters for PHA biosynthesis and a comparative PHA genetic organization for Bacillus sp. BNPI-92 and other closely related strains from RAST server annotation for a and b using online tools that are available at http://rast.theseed.org/FIG/rast.cgi. A deep red colored arrow indicates a phaC gene that is supposed to encode protein used for PHA polymerization. a is a phaA gene or nucleotide sequences (744 bp) encoded by acetyl-CoA thiolase and b is a phaC gene or nucleotide sequence (1086 bp) encoded by polyhydroxyalkanoic acid synthase. It is a conserved nucleotide sequence. And phaB gene (arrow with number 2) encodes acetoacetyl-CoA reductase. phaB and phaC genes are located on the same operon. Note that phaB gene sequence is not shown. c Three-dimensional structure of phaC protein with 319 residues (residue of protein structure for phaC gene) and presumed for PHA biosynthesis in Bacillus sp. BNPI-92. It has 319 residues. Its 88% has been modeled with 100.0% confidence by the single highest scoring template. The given protein residue resembles class I polyhydroxybutyrate synthase that was derived from Cupriavidus necator in terms of its structure. The three-dimensional structure of protein was predicted using PHYRE2 which is an online tool that is available at http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index (d). Evolutionary relationships for Bacillus sp. BNPI-92

Fig. 8

PHA biosynthesis promoting genetic organizations (gene loci) in Bacillus sp. BNPI-92 and other closely related type strains genome ( Bacillus cereus ATCC 14579 T (AE016877), Bacillus paranthracis Mn5 T (MACE01000012), Bacillus thuringiensis ATCC 10792 T (ACNF01000156), and Bacillus antrics Amen T (AE016879). These strains harbored PHA biosynthesis gene with nearly similar nucleotide sequence size (nts) such as croR (predicted for PHA metabolisms), PhaQ (a PHB-responsive repressor controlling expression of phaP ), phaP (predicted for phasin biosynthesis), phaA (suggested for acetyle fermentation), phaB (predicted for butyrate fermentation), phaR (suggested as class IV PHA synthase or polyhydroxyalkanoate synthesis repressor phaR ), and phaC (predicted as PHA polymerization). phaA was located on a separate gene and suggested for encoding 3-ketoacyl-CoA thiolase or acetyl-CoA acetyltransferase protein. phaJ [5] (encoding enoyl-CoA hydratase) and acsA gene [49] (encoding acetoacetyl-CoA synthetase) predicted for PHA biosynthesis were located on the same operon. Isolog butyryl-CoA dehydrogenase enzyme is a gene which is unidentified

Fig. 9

Metabolic pathway for PHA-producing Bacillus cereus BNPI-92 strain obtained from an area of plastic waste accumulation. A letter in green color box (a–v) indicates the most likely enzymes involved in metabolic activities. The letter in violet color box signifies possible gene involved in coding enzyme. List of enzymes: a not predicted in this study, b formate C-acetyltransferase [EC:2.3.1.54] (transferases), c acetyl-CoA C-acetyltransferase [EC:2.3.1.9] (transferases), d acetoacetyl-CoA reductase [EC:1.1.1.36] (oxidoreductases), e polyhydroxyalkanoate synthase subunit phaC [EC:2.3.1.-] (transferases), f butane monooxygenase alpha subunit [EC:1.14.13.230] (oxidoreductases), g butanol dehydrogenase [EC:1.1.1.-] (oxidoreductases), h acetaldehyde dehydrogenase (acetylating) [EC:1.2.1.10] (oxidoreductases), i butyryl-CoA dehydrogenase [EC:1.3.8.1] (oxidoreductases), k 3-hydroxybutyryl-CoA dehydratase/vinylacetyl-CoA-Delta-isomerase [EC:4.2.1.1205.3.3.3] (lyases), l butyrate kinase [EC:2.7.2.7] (transferases), m phosphate butyryltransferase [EC:2.3.1.19] (transferases), n 2-hydroxyglutarate dehydrogenase [EC:1.1.99.2], p glutaconate CoA-transferase, subunit A [EC:2.8.3.12] (transferases), q unidentified, r glutaconyl-CoA decarboxylase subunit alpha [EC:7.2.4.5], s 3-hydroxybutyryl-CoA dehydratase [EC:4.2.1.55] (lyases), t 3-hydroxybutyrate dehydrogenase [EC:1.1.1.30] (oxidoreductases), u hydroxybutyrate-dimer hydrolase [EC:3.1.1.22] (hydrolases), and v poly(3-hydroxybutyrate) depolymerase [EC:3.1.1.75] (hydrolases). List of few predicted gene in these pathways: a unidentified, b pf1D or it could be pflB [50] gene , c atoBgene, and d phbB gene. However, others were designated as phaA [–53], orphbA [54, 55] gene, d its phaB [56, 57] or phbB [58] gene, e predicted as phaC or phbB gene , f unidentified in pathway, g it could be yugJ [59], h unidentified in the KEGG-KASS pathway, i gene encoding isologs of butyryl-CoA dehydrogenase enzyme are unidentified. It could be Swol_1933 [60], k croR gene, l unidentified in pathway, n unidentified gene, p it could be gctA , r it could be gcdA (s) croR gene, t it could be Bdh2 [61] gene, u unidentified, and v phaZ gene. These metabolic pathway and respective KEEG number database were collected from KEGG [62]