Whole Genome Sequence of biosurfactant producing Bacillus tequilensis

2.1 Identification of Biosurfactant producing Bacillus tequilensis

Majority of biosurfactants are produced by the microbes such as Pseudomonas genus followed by Bacillus and Acinetobacter respectively[11]. In a previous investigation, a novel strain of biosurfactant producing Bacillus tequilensis strain ANSKLAB04 [12] was identified using 16S rRNA gene sequencing by Sanger dideoxy sequencing method followed by the phylogenetic assessment. The strain was isolated from Chilika lake, a brackish water lagoon, spread over the Puri, Khurda and Ganjam districts of Odisha state on the east coast of India[12]. By conducting several biochemical tests such as Haemolysis test, oil spreading test, CTAB agar plate test and Drop collapse test, we concluded that Bacillus tequilensis produces biosurfactants [12] and the novel isolate was deposited in Genbankwith Accession number KU529483.

2.2 Bioanalyzer profile

The DNA isolation was performed using Phenol/Chloroform(PCl) genomic DNA extraction method[12].The bioanalyzer profile of the prepared WGS library showed fragments in a size range of 300–600 bp. The effective insert size of the library is 180–480 bp flanked by adaptors having combined size of ~ 120 bp. Based on the fragment distribution and concentration, the library was suitable for sequencing on Illumina platform [Fig. 1].

2.3 Genome Representation

The complete genome of Bacillus tequilensis consists of a single circular chromosome of 4,478,749 bp with an average G + C content of 46.33% (Table 1 and Supplementary Table 1). The 4492 predicted coding ORFs covers 87% of the complete genome, and each ORF has a moderate length of 283 aa(Supplementary Table 2). Among these, 1,347, i.e. 67.4% assigned as putative functions, 258, i.e. 12.9% matched to sustain hypothetical coding sequences of an anonymous function, and the rest 394, i.e. 19.7% shows no similarities to known genes[Table 2].

The variations in nucleotide frequencies across the whole genome sequence was investigated using a non-overlapping active platform and by framing three indices of nucleotide frequency: G + C%, (G + C)/(A + T + C + G), divergence from [A]=[T], (A-T)/(A + T), and divergence from [C)=(G), (C-G)/(C + G). These 3 indices are, by represent, pairwise-independent and summarize relative nucleotide frequencies without loss of information. Because of their very low frequency, ambiguous nucleotide bases were not taken into account. The SD (standard deviation) for the 3 indices are given by

Normal distribution approximation was used as the total number of bases were large. The strand analyzed here was the 5' to 3' strand clockwise on the genetic map. A window size of 1 kb was used.From the inside: green and red bars represent RNA sequences on positive and negative strand respectively. Circle 1, represents G + C content (window size: 10Kb) higher and lower than 45%, where red represents higher and green represents lower. Circles 2: - represents GC skewness, where green and red represents positive and negative value respectively [Fig. 2].

2.4 Gene ontology and biological annotation

The gene ontology analysis concluded that 18.99% of genes in Bacillus tequilensisbelonged to transferase activity, 13.55% of genes belonged to kinase activity, 9.3% of genes were involved in ATP binding, 9.3% genes were involved with hydrolase activity, 6.91% genes were involved in methyltransferase activity, 5.98% of genes were associated with lipase activity, 5.98% of genes were involved in oxidoreductase activity, 4.9% of genes were in lyase activity, 3.05% genes were involved in peptidase activity, whereas only 2.79% genes were involved in cell division, 2.9% were in carbohydrate transport, 7.7% were in ribose production, and only 3.8% genes were involved in viral capsid [Fig. 3].

2.5 Subsystem Classification

Genes obtained from the whole genome of Bacillus tequilensishas been used for the classification of the subsystem. Subsystems were categorized based on the cofactors, cell wall, virulence metabolism, potassium metabolism, membrane transport, iron acquisition and metabolism, RNA metabolism, cell division and cell cycle, motility and chemo taxis, fatty acids, lipids and isoprenoids, nitrogen metabolism, etc were discussed in [Supplementary Table 3 ][Fig. 4].

2.6 Metabolism of biosurfactant producing genes

Bacillus tequilensisproduces a biosurfactant that belongs to the class of lipopeptides having excellent emulsifying properties and were capable of reducing the surface tension of water to a significantly lower value. The genes associated with producing biosurfactant are listed in [Table 3]. Among the several different classes of Biosurfactant producing bacteria genera, the members of the genera Bacillus or Pseudomonas, due to their wide range of applications and resourcefulness can be more often is used. Bacillus species are phenotypically and genotypically heterogeneous. Based on several investigations, a unique inhabitant of Bacillus sp. found at the marine site such as B. subtilis, B. licheniformis, B. cereus, B. amyloliquefaciens, B. pumilus, and B. mycoides. Bacillus subtilis produces lipopeptide biosurfactant called surfactin, which is coded by four ORFs named as SrfA, SrfB( also known as ComA), SrfC, and SrfD. The sfp gene considered as an essential component of peptide synthesis systems and plays a major role in the regulation of surfactinbiosynthesis and gene expression. Srf gene amplification is at 268 bp whereas the expression of the sfp gene amplified at 675 bp [13]. The peptide synthesise for an amino acid moiety of surfactin is encoded by four ORFs in the srfA operon namely SrfAA, SrfAB, SrfAC and SrfAD /SrfA-TE and also contains comS gene lying within the out-of-frame with the srfB [14]. PorobS et al 2013 and Nakano Met al., 1992 isolated SrfA gene from Bacillus amplified at 580 bp and authors concluded the biological significance of SrfA gene in biosurfactant production [15–17]. From Bacillus tequilensis we identified the SrfA which involved in biosurfactant production and the sequence of the SrfA(242 aa) was deposited in GenBank with accession MUG02427.1.

Besides, lichenysin is another lipopeptide biosurfactant produced by B. licheniformis coded by lichenysin operon (LchA) and comprises of four peptide synthetase genes: LicAA, LicAB, LicAC, and LicAD. In another study, the authors isolated genesfp (Phosphopantetheinyl transferase 224 amino acids) and mapped at 4 kb downstream to operon srfAand the authors also concluded it is essential for the post-translational changes to surfactin synthetase in microbes [15–16]. In this study we have identified sfp gene from Bacillus tequilensis and the sequence of the sfp (Phosphopantetheinyl transferase 224 amino acids) was deposited in GenBank with accession MUG02422.1.

Moreover, two operons, srfA and pps were found to be present in UMX-103 and B. subtilis 168 strains only involved in biosurfactant synthesis. The srfA operon contains four genes such as srfAA, srfAB, srfAC, and srfAD and the operon pps contains four genes named as ppsB, ppsC, ppsD, and ppsE. The genes, rmlA, rmlB, rmlC and rmlD are only present in UMX − 103 strains whereas, sigA, DnaK and LytR are present specifically in Bacillus strain. Besides, the genes comA, comP, rpoN, abrB and ResD are presented in both UMX-103 and B. subtilis 168 [18]. Based on the above literature biological annotation, we have identified DnaKandLytRgenes from Bacillus tequilensis and the sequence were deposited in NCBI with accession MUF99480.1 and MUG01692.1 respectively.

Pseudomonas species required Plasmid-encoded- rhlA, B, R and I genes of rhl quorum-sensing system for production of glycolipid biosurfactants as well as also involved in the production of rhamnolipids in a heterologous host. Iturin A is an antifungal lipopeptide biosurfactant produced by certain Bacillus subtilis strains such as B. subtilis RB14 is composed of four ORF namely ituD, ituA, ituB, and ituC, whose disruption leads to specific deficiency in iturin A production. The three genes of arthrofactin operon of Pseudomonas namely arfA, arfB and arfC encode ArfA, ArfB and ArfC containing two, four and five functional modules respectively required for condensation, adenylation and thiolation. Besides, Amphisin is produced by Pseudomonas sp. DSS73 require gacS and amsY genes for the production of biosurfactant as these genes are mutants defective in the genes. Amphisin synthesis is regulated by the gacS gene as the gacS mutant regains the property of surface motility upon the introduction of a plasmid. Moreover, genes dnaK, dnaJ and grpE positively regulate the biosynthesis of putisolvin [14].Putisolvin biosynthesis genes such as dnaK, dnaJ and grpEfrom Bacillus tequilensis were identified and the sequence were deposited in Genbank with accession MUF99480.1 MUF99481.1, MUF99479.1 respectively.

Acinetobacter species produces high molecular weight biosurfactants - Emulsan and Alasan with the involvement of gene. AlnA, AlnB and AlnC are essential for Alasan biosynthesis whereas wza, wzb, wzc, wzx, and wzy is required for Emulsan biosynthesis. For the production of fungal biosurfactant, emt1 and cyp1 are the two genes involved in the synthesis of these glycolipids and fb1 and hfb2 genes regulating the synthesis of hydrophobin [14]. Thus, gene plays a major role in the biosynthesis of various microbial surfactants, and hence the role of molecular genetics and gene regulation mechanisms in the production of biosurfactant is essential. In this study, we have identified biosurfactant producing genes and corresponding orfs of Bacillus tequilensis such gene SrfAD, SrfAC, SrfAA and the sequence of the same was deposited in Genbank database with accession MUG02427.1, MUG02428.1, MUG02429.1,MUG03515.1 respectively.

2.7 Biosurfactant / Lipopeptide metabolism of Bacillus tequilensis

Considering the biosurfactant producing genes described in various literatures, we classified the genes of Bacillus tequilensis based on the established efficient biosurfactant activity and broad applications. Biosurfactant is proven to be promising; possessing unique properties of low toxicity and higher biodegradability. In the present investigation we constructed a pathway which describe the biosurfactant metabolism of Bacillus tequilensis[Fig. 5]. The lipopeptide synthesised constitutes of a long chain of fatty acid along with glutamate acid (Glu), leucine (Leu), aspartic acid (Asp) and valine (Val). The synthesis is non-ribosomal by a large multienzyme peptide, non-ribosome peptide synthases (NRPS). The peptide synthetase required for an amino acid moiety of surfactin is encoded by four open reading frames in the srfA operon namely SrfAA, SrfAB, SrfAC and SrfAD or SrfA-TE. SrfA, SrfB, SrfC and srfD constitute the four main enzymes for surfactin formation. SrfD is the most important enzyme as it initiates the surfactin formation. Sfp gene is 224 amino acid long- present downstream to srfA operon- also plays an important role as it is required for the posttranslational modifications to surfactinsynthetase.

Different modules have been marked based on different pathways involved for the synthesis of biosurfactant, such as glycolysis, TCA cycle, NADPh generation, amino acid biosysnthesis, fatty acid synthesis and synthesis of surfactin. Seven modules represent the different pathways required for production of glutamate acid (Glu),leucine (Leu), aspartic acid (Asp) and valine (Val). The precursors for biosynthesis of Val/Leu, Glu/ Asp, and fatty acids are the product of glycolysis and TCA cycle such aspyruvate, 2-oxo-glutarate, oxaloacetate, and acetyl-CoA. The genesof Bacillus tequilensisinvolved in the utilization of sucrose,including sacP, murPand sacA, which encode a sugar transporter, permease, and sucrose-6-phosphate hydrolase,were identified and the sequence was deposited in genbank with accession MUF99868.1,MUG00557.1, MUG01465.1 respectively. The NADPH generation and pentose are produced by pentose phosphate pathway catalysed by zwfand GNDA enzyme.

The biosynthesis of Glu, Asp, Val, and Leu, are considered as the intrinsic components of surfactin.Glu/Asp are synthesised by aspartate aminotransferase such as AspB and YhdR were identified from Bacillus tequilensis and the sequence were deposited to genbank using accession MUF99794.1 and MUF99877.1respectively. The efficient fatty acid biosynthesis pathway determines efficient surfactin production. The building precursor acetyl-CoA initiates the biosynthesis of fatty acid. The biosynthesis of surfactin is catalysed through NRPS, initiated from the condensation of fatty acids andGlu. Other constituent amino acids are assembled through the NRPS multi-enzyme complex, comprising adenylation,condensation, and thiolation domains responsible for the activation of amino acids and peptide chainelongation.

2.8 Genome evolution of B. tequilensis

The enormous genomic data obtained from sequencing of Bacillus tequilensisANSKLAB04 was aligned against existing top 20 homologous species of bacillus in the NCBI database. The comparison of the number of unique genes and common genes were analysed with top 6 homologus species of bacillus such asB. subtilis, B. vallismortis, B. tequilensis (KCTC 13622), B. halotoleran and B. mojavensis[Table 4][Fig. 6].

The phylogenetic tree was constructed from the top 20 homologous bacillus species obtained from a blast search. The orange colors show the gene family expansion and grey color indicates the gene family contractions between the bacillus species. The corresponding proportions among the total changes are shown in the same colors in the pie chart. Implied divergence dates (in millions of years) are indicated at each node in blue. The most recent common ancestor (MRCA) and the blue color indicates the conserved gene family among the various species of Bacillus[Fig. 7].

2.9 SSR Identification

Microsatellites or simple sequence repeats (SSRs) also known as short tandem repeats (STRs) are broadly used as PCR-based markers which can be helpful for characterization of population genetics, genetic mapping, phylogenetics, genome mapping, population genotyping, quantitative train loci analysis, genome comparisons, and markers assisted breeding etc.The present investigation predicted the accurate repeated SSRs with repeat unit length between 1 and 10 bp across the whole genome (528 scaffolds. i.e. 4478749 bp) of the whole genome of Bacillus tequilensisANSKLAB04. K-mer distribution analysis found a total of 25 motif 2-mer, i.e. 86.2%, 3 motifs 3-mer i.e. 10.34% and motif 5 -meri.e. 3.4% in the whole genome of Bacillus tequilensisANSKLAB04 [Fig. 8] [Fig. 9].

SSR composition was studied from each scaffold. The present investigation classified the SSR based on monomer, dimer, trimer, tetramer, pentamer, hexamer as P1, P2, P3, P4, P5, P6 respectively. The composition and repetition number of P1 to P6 lengthwise were studied and provided in [Table 5–6]. A total of 32 SSR were found from the 528 sequences of Bacillus tequilensisANSKLAB04. 27 Sequences were identified with SSRs [Table 7].A total of 27 P1 were found which recorded as the highest and P6 found the lowest [Fig. 9].

2.10 SNP and Indel Discovery

Our Indel discovery strategy involved mining insertion and deletion polymorphisms from DNA sequencing traces that originally were generated by genome centres for SNP discovery. The obtained mass sequenced data of Bacillus tequilensisANSKLAB04 were used to search for genetic variation against existing homologous biosurfactant producing bacteria from NCBI. The present investigation used the existing to 5 homologous genomes of bacteria such as Bacillus tequilensis KCTC 13622, Bacillus subtilis, Bacillus mojavensis, Bacillus vallismortis, Bacillus halotolerans. The number of mapped sites per sample, mapping coverage, the total number of reads, number of mapped reads, overall mapping ratio, number of mapped bases, and the average alignment depth was calculated.Table 8represents the statistics of Bacillus tequilensison comparison with 5 existing homologous with reference bacterial genome which includes Bacillus tequilensis(KCTC 13622), Bacillus halotoleran, Bacillus subtilis, Bacillus mojavensis and Bacillus vallismortis. The number of total reads in all the reference genome is 6,229,938 which are constant in all reference bacteria. The mean depth indicates the number of reads, on average, are likely to be aligned at a given reference base position on comparison with Bacillus tequilensis. However, Bacillus subtilis is having 90.39% of mapped read, 786,017,247 mapped bases and 186.45 mean depth which is highest among the other indicative of better analogy and susceptibility. On the other hand, Bacillus mojavensis with reference length 3,957,021, mapped reads 73.19%, mapped bases 555,891,216 and mean depth 140.48 showing the least compatibility with the Bacillus tequilensis.

After removing duplicates with Sambamba and identifying variants with SAMTools, information of each variant were gathered and classified by chromosomes or scaffolds. Table 9 shows the summary of the variant calling of Bacillus tequilensis ANSKLAB04 against other existing genomes in the database.

Table 9 represents the summary of Variant Calling of Bacillus tequilensis against existing top 5 homologous references bacterial genome which includes Bacillus tequilensis (KCTC 13622), Bacillus halotoleran, Bacillus subtilis, Bacillus mojavensis and Bacillus vallismortis. Comparison of the whole-genome sequence of Bacillus tequilensis on comparison with a reference reveals the number of markers that include single nucleotide polymorphisms (SNPs), inserted and deleted sequences. Figure 10represents the graphical representation of SNPs and INDEL in which Bacillus halotolerans is having the highest number of SNPs i.e. 347,175 [Figure 10D] whereas Bacillus tequilensis KCTC 13622 is having the maximum number of insertions and deletions i.e. 841 and 653 respectively [Figure 10A]. Meanwhile, Bacillus subtilis were having less number of SNPs, insertions and deletions [Figure 10B].

4.1 Sample Collection and DNA Isolation

Water samples were collected from oil-contaminated sites of Chilika lake, Odisha, India (latitude and longitude: 19.8450 N 85.4788 E), a largest brackish water lagoon in India. Various organisms were isolated and purified on culture plates and were then enriched in the mineral salt medium (MSM). MSM gives the nutrient condition for the production of biosurfactants by the organisms which were then screened for their biosurfactant production by various screening tests and the emulsification index was calculated. Identification of organisms was performed based on biochemical, macroscopic and microscopic characters. The organism with the best emulsification index was then subjected to optimization for the production of biosurfactant for the factors affecting the production. Optimization was studied with the emulsification index calculated with each affecting factor. In a previous study, this organism was then subjected to 16S rRNA sequencing for the identification of the genus and species [12]. The DNA was isolated by Phenol/Chloroform (PCl) genomic DNA extraction method [12][19]. The bacterial cell pellet obtained after centrifugation was subjected to DNA isolation. The DNA concentration and purity were checked with nanodrop spectrophotometer and qubit fluorometer. ImageQuant software was used to analyse the gel image documentation [20].The sample was found to have suboptimal concentration and gave an intact band when running on gel QC against Hind III digested lambda ladder. The gel image demonstrates the intact Bacillus tequilensis. Whole genome having a concentration of 54.9 ng/ul. L represents Hind III Digested lambda ladder. '1 - SO_4915_Bt1' represents the Bacillus tequilensis sample code. [Fig. 13].

4.2 Materials used in the study

Whole Genome Sequencing kits such as NEXTFlex DNA Sequencing Kit (Cat # 5140-02), NEXTFlex DNA Barcodes – 48 (Cat # 514104), HighPrep™ PCR (Magbio, #AC-60050),High Sensitivity Bioanalyzer Chips (Agilent, #5067 − 4626), Nuclease free water (Ambion, #AM9939), Covaris™ S220 System (Life Technologies, #4465653), Covaris™ microTUBE AFA (Life Technologies, #520045), Low Melting Agarose (Invitrogen, #16520100), MinElute Gel Extraction kit (QIAGEN, #28604), 50X TAE Buffer (MP Biomedical, Cat #TAE50 × 01), Qubit® dsDNA HS Assay Kit (Invitrogen, Cat # Q32854) were used in the present investigation.

4.3 Library Preparation and Genome Sequencing

Library preparation was performed using NEXTFlex DNA library protocol outlined in “NEXTFlex” DNA sample preparation guide (Cat # 5140-02). In brief, genomic DNA was sheared to generate fragments of approximately 300-500bpin a CovarismicroTube with the E220 system (Covaris, Inc., Woburn, MA, USA). The fragment size distribution was checked using Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA) with High Sensitivity DNA Kit (Agilent Technologies) according to the manufacturer’s instructions. The resulting fragmented DNA was cleaned up using HighPrep beads (MagBio Genomics, Inc, Gaithersburg, Maryland). These fragments were subjected to end-repair, A-tailing, and ligation of the Illumina multiplexing adaptors using the NEXTFlex DNA Sequencing kit as per the manufacturer’s instruction [21].

The resulting ligated DNA was cleaned up using HighPrep beads (MagBio Genomics, Inc, Gaithersburg, Maryland)and size selected (400–600 bp) on 2% low melting agarose gel and cleaned using MinElute column (QIAGEN, India). These adapter-ligated fragments were subjected to 10 rounds of PCR (denaturation at 98˚C for 2 min, cycling (98˚C for the 30S, 65˚C for 30S and 72˚C for 1 min) and a final extension at 72˚C for 5 min) using primers provided in the NEXTFlex DNA Sequencing kit(Perkin Elmer). The PCR products were purified using HighPrep beads. Quantification and size distribution of the prepared library was determined using Qubitflourometer (Table 18) and the Agilent High Sensitivity DNA Kit (Agilent Technologies) respectively according to the manufacturer’s instructions (Fig. 14). Illumina Paired-end sequencing was performed using NextSeq 500: 150*2. The following adapters were used for sequencing (Illumina, Inc)[21].

Adapter details: Universal Adapter 5’ AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT 3’

Adapter, Index 5’

GATCGGAAGAGCACACGTCTGAACTCCAGTCAC[INDEX]ATCTCGTATGCCGTCTTCTGCTTG 3’

4.3 Whole Genome De-novo Assembly and Analysis

The obtained sequence raw reads were checked for quality control using FASTQC tool [22]. The quality of the raw reads was checked through the various modules provided by the FASTQC tool. Among the modules, per base sequence quality and tile sequence quality modules were studies to validate the quality of the data for further analysis. The low-quality reads were excluded from the analysis using Trimmomatic (v0.36)[23]. The filtered De-novo assembly of Illumina paired-end data was assembled using SPAdes - v3.13.0 genome assembler - an open-source algorithm for De-novo assembly [24]. SPAdes assembler is intended for de-novo assembly after error-correction of sequenced reads. Assembled contigs were further scaffolded using SSPACE program [25]. Genome map was constructed using Circos [26].

4.4 Whole Genome annotation and GO analysis

NCBI Prokaryotic Genome Annotation Pipeline (PGAP) version 4.8 was used to annotate the whole genome sequence of [27]. Pathway Analysis was done by using KAAS Server. Bacillus subtilis subsp. Subtilis 168 was taken as a reference organism for pathway analysis using KAAS server [28]. The functions of the predicted ORFs were categorized by comparison with the COG database [29]. Simple Sequence Repeats (SSR) was identified in each transcript sequence using MISA Perl script [30].

4.5 Variant Calling and Variant annotation

Variant Calling of Bacillus tequilensiswas performed by aligning with the top 5 existing homologous reference bacterial genome which includes Bacillus tequilensis(KCTC 13622), Bacillus halotoleran, Bacillus subtilis, Bacillus mojavensis and Bacillus vallismortis. The present investigation used BWA(Burrows-Wheeler Aligner)-MEM for alignment of Bacillus tequilensisagainst top 5 homologous genomes[31]. During mapping, duplicated reads can falsely cause erroneous data to stand out. To prevent such error, Sambamba tool was used to remove the duplicate reads [31]. Duplicate reads are identified using mapping information such as start position, and CIGAR string [32]. SAMTools was used to manipulate the SAM/BAM files that come out as a result of mapping [33]. In resequencing analysis, it is especially used for finding out variant information by calculating genotype likelihood from every position within the sample of analysis. Variant annotation was performed using SnpEff (v4.3t)[34]. SnpEff annotates the possible effects (on genes) that can be caused by variants identified through mapping.The present study used SnpEff to generate the Genes and transcripts affected by the variant, Location of the variants and the information on how the variant affects the protein synthesis (e.g. generating a stop codon).