Introduction

The gut microbiota is a complex microecosystem that harbours around 100 trillion microbes derived from over 2,000 diverse species [1]. It plays a vital role in maintaining the energy balance, regulating gut pH and developing gut cells in the host body [2]. Some species of the gut microbiome including Bacillus, Lactobacillus, Lactococcus, Bifidobacterium, Streptococcus, Enterococcus, and Saccharomyces cerevisiae have been classified as probiotics [2]. Probiotics are microorganisms that provide a health benefit to the host when administered at specific concentrations [3]. They have also been reported to produce antimicrobial compounds having bio-therapeutic potential [4] and are involved in the immunological development and intestinal barrier improvement in the host [1].

In order to be recognized and used as a probiotic, a microbe must possess several characteristics such as adherence, a range of activities against pathogenic microbes, tolerance to harsh conditions, etc. to survive the upper intestinal tract and reach the site of action. The Bacillus spp., due to their spore-forming ability and higher rate of secondary metabolism are promising probiotic candidates. The endospore structure is able to survive harsh GIT environment, such as high temperatures, low pH and high bile salts, and the range of active substances synthesized confer anti-cancer and antioxidant properties to the strain [2]. These characteristics are responsible for the higher success rate of Bacillus in colonizing the GIT as compared to other genera [5]. Bacillus subtilis, first described in the 19th century, is a Gram positive, aerobic, fast growing, spore forming bacterium [6]. Bacillus subtilis strain 168 was one of the first bacteria whose genome was fully sequenced and remains as one of the best annotated genomes [7,8]. B. subtilis as a probiotic offers multiple health benefits such as anti-allergic effects, increased immunity and reduced bone loss in postmenopausal women [9]. Moreover, Bacillus isolates also produce a wide range of antimicrobial compounds, including lipopeptides and Bacteriocin-Like Inhibitory Substances (BLIS) [10]. Several species of the genera including B. subtilis, B. clausii, B. licheniformis, B. coagulans, B. polyfermenticus, B. cereus and B. pumilus are widely used as probiotics [6].

Even though consumed regularly, there are a few risks associated with some B. subtilis strains, such as the presence of genes responsible for enterotoxin and biogenic amine synthesis, possibility of AMR gene transfer, and cytotoxicity. Therefore, before a strain could be used as a probiotic, it is essential to evaluate its safety and efficacy by assessing parameters like presence of potential virulent or pathogenic factors, toxin and biogenic amine production, and antibiotic resistance. Wholegenome sequencing (WGS) technology has been employed not only in the characterization of promising probiotic strains but also for the evaluation of its safety aspects in an efficient manner. [1116]. Bacterial strains which differ by even a single nucleotide can be identified by WGS. Particularly, Food and Drug Administration (FDA) has started emphasizing on WGS to identify food borne pathogens and prevent illness. EFSA (European Food Safety Authority) also recommends WGS data for the approval of microbes as feed or food additives [11].

In this study, the WGS of B. subtilis strain PLSSC was performed using Illumina and Nanopore technology. The strain was identified and confirmed by analysing the 16S rRNA gene sequence, mol% G+C content, and average nucleotide identity (ANI). Genome analysis suggested that B. subtilis strain PLSSC is safe and has genes essential for the probiotic characteristics.

Materials and Methods

Extraction of genomic DNA and purification

GeneAll Kit (GeneAll, Seoul, South Korea) was used for the isolation and purification of B. subtilis PLSSC genomic DNA. Bacterial culture (2 ml) was grown overnight and centrifuged at 8000×g for 2 minutes. After centrifugation, pellet was collected and re-suspended in 180 µl of GP buffer (containing 30 mg ml−1 of lysozyme) and subsequently incubated at 37 °C for 30 minutes. After incubation, RNAase and Proteinase K were added to the suspension. Prior to transfer to a binding column, 200 μl of absolute ethanol was added to the lysate. Column was washed using a wash buffer and afterwards, elution of DNA was carried out with 100–200 µl of sterile Milli-Q water. After isolation, genomic DNA was assessed in terms of quantity and quality by NanoDrop-2000 (ThermoFisher Scientific, USA), Qubit, and agarose gel electrophoresis.

Whole genome sequencing and assembly

For the preparation of Illumina and Nanopore Whole Genome Sequencing (WGS) libraries, purified genomic DNA from B. subtilis strain PLSSC was used. Illumina WGS library was prepared by using Illumina-compatible SureSelect^QXT (Fig. 1) whole genome library prep kit (Agilent, Santa Clara, CA, U.S.A.).

Figure 1: Workflow for illumina SureSelect^QXT library preparation.

25 ng DNA was subjected to fragmentation and adapter tagging using Sure Select QXT enzyme. HighPrep PCR beads were used to purify fragmented and adapter-tagged DNA. Amplification and indexing of fragmented and adapter-tagged DNA were carried out using 6-cycles of PCR. Again, HighPrep beads (MAGBIO, MD, USA) were used to purify PCR products. Purified PCR products were subjected to library quality control check. Quantification of Illumina-compatible sequencing library was done by Qubit fluorometer (Thermo Fisher Scientific, MA, USA) and fragment size distribution was analysed on Agilent TapeStation (Table 1 and Fig. S1).

The tape station profile of Illumina library reveals that the fragment size ranges from 187 to 1149 bp. However, larger proportion of the Illumina-compatible sequencing library had a fragment size ranging from 200 to 700 bp. Taking into account the combined adapter size i.e. approximately 120 bp, the effective user-defined insert size was 80 to 580 bp, which was obtained with optimal concentration. This ensured that the library was suitable for Illumina sequencing to get the desired amount of sequencing data. A total of 2,490,092 (R1 + R2) reads were generated for B. subtilis strain PLSSC on Illumina MiSeq platform.

1 µg of genomic DNA from B. subtilis strain PLSSC was used for the generation of long read library for Nanopore sequencing. Firstly, the genomic DNA was end-repaired with NEBnext ultra II end repair kit (New England Biolabs, MA, USA), followed by

a clean-up with 1x AmPure beads (Beckmann Coulter, USA). Library preparation was then performed using a native barcoding kit (NBD103) wherein barcodes were ligated using NEB blunt/TA ligase (New England Biolabs, MA, USA), followed by a clean-up with 0.5x AmPure beads (Figure 2).

Figure 2: Native barcoding library preparation.

Nanopore sequencing was executed on GridION X5 (Oxford Nanopore Technologies, Oxford, UK) using SpotON flow cell R9.4 (FLOMIN106) in a 48 hour sequencing protocol on MinKNOW 2.1 v18.05.5. Albacore v2.3.1 was used for base-calling Nanopore raw reads in ‘fast5’ format to ‘fastq’ format.

MaSuRCA Hybrid Assembler [17] was used for generating a hybrid assembly of Illumina and nanopore reads. Gene annotation of assembled genome was carried out by NCBI Prokaryotic genome Annotation Pipeline [18].

Table 1: Description of the library.

Calculation of average nucleotide identity

Average Nucleotide Identity (ANI) is used to measure the relatedness between two genomes. ANI calculation of the assembled genome of B. subtilis strain PLSSC against the genome of B. subtilis KCTC 3135 was performed according to Goris et al. [19] and calculated using ANI calculator (http://enve-omics.ce.gatech.edu/ani/) with default parameters.

Identification of antibiotic resistance and virulence factor genes

For the identification of antibiotic resistance associated genes, a homology-based search was carried out between the assembled genome of B. subtilis strain PLSSC and Comprehensive Antibiotic Resistance Database (CARD) [20]. For the identification of hits, BLASTX was used with similarity >30%, coverage >70% and e-value <1e-02. Additionally, the assembled genome was also compared with the COG database [21] to identify gene function.

Virulence Factor Database (VFDB) [22] having 3072 sequences in the core database was used to search for the virulence factors genes. Search was performed using BLASTX and only those hits which showed a similarity >30%, coverage >70% and e-value <1e-02 were taken into consideration.

Identification of biogenic amine producing genes

Genes involved in biogenic amine production, mainly amino acid decarboxylases were searched in the assembled genome of B. subtilis strain PLSSC as described by Salvetti et al. [23]. Protein sequences of the short-listed biogenic amine producing genes (amino acid decarboxylases) were downloaded from the Uniprot database and BLASTX was performed between the assembled genome and biogenic amine producing proteins.

Assessment of Genomic Stability

Stability of the assembled genome was assessed in accordance with Salvetti et al. [23]. Presence of insertion sequences (IS), prophage sequences and clustered regularly interspaced short palindromic repeats (CRISPR) sequences in the genome was investigated. PHASTER, a web-based server was used to identify prophage sequences in the assembled bacterial genome [24]. Mobile elements were searched using ISfinder (web-based software) and ACLAME database (version 0.4). For screening CRISPR sequences, CRISPRCasFinder was used [25].

In silico mining of probiotic genes

The assembled genome of B. subtilis strain PLSSC was assessed for the genomic features which contribute to the probiotic properties such as adhesion to gut mucosa, acid tolerance, bile salt tolerance and environmental stress resistance [26]. All the predicted protein sequences (n=4296) were annotated by submitting them to Batch-CD Search web service available in the conserved domain database. In the Batch-CD search, the Pfam database containing 19178 position specific scoring matrices (PSMMs) was selected for functional annotation of predicted proteins.

Results

Assembly of B. subtilis strain PLSSC genome and its features

De novo WGS of B. subtilis strain PLSSC yielded a single scaffold of 4,204,670 bp in size. Annotation of the assembled genomic sequence by the NCBI Prokaryotic Genome annotation pipeline resulted in 4296 coding sequences (CDS), 86 tRNAs, 30 rRNAs and 5 ncRNA, which were in accordance with the annotation described for B. subtilis KCTC 3135 (https://www.ncbi.nlm.nih. gov/nuccore/CP015375.1/). The predicted protein coding genes were annotated against all the proteins (228,533) listed in the Uniprot database for bacteria (Uniprot Consortium, 2014). We observed that, out of 4,296 CDS from the assembled genome, 4,256 got significant hits (identity >30% and e-value <=1e-05) against the Uniprot Bacterial protein database. For these 4,256 proteins, a gene ontology classification (molecular function, cellular component and biological process) was also carried out and represented in the form of a pie chart (Figure 3). Gene ontology analysis indicates that 39.25%, 45.49% and 15.26% of the proteins from the assembled genome represented for molecular function, cellular component and biological process respectively. Major pathway groups associated with the predicted proteins and their functions are indicated in Figure 4 and 5.

Figure 3: Gene ontology (GO) association of predicted protein-coding genes.

Using the same criteria, predicted protein coding genes from the assembled genome were also annotated with Clusters of Orthologous Groups (COG) database (http://www.ncbi.nlm.nih.gov/COG/) [21]. It was observed that, out of 4,296 CDS from the assembled genome, homology was found against the COG database for 4,091.

Figure 4: Pathway abundance of predicted proteins. Only top 30 pathways have been shown.

Figure 5: Pathway function associated with predicted proteins.

Taxonomic analysis

As mentioned before, the assembled PLSSC genome contains 30 rRNA genes. Further classification of these 30 rRNA genes revealed that there were equal number (10) of 5S rRNA, 16S rRNA and 23S rRNA genes. Upon alignment to the SILVA 16S Database [27], ten 16S rRNA sequences exhibited significant similarity among themselves. Moreover, phylogenetic tree constructed from 16S rRNAs showed that the assembled strain PLSSC is closely related to Bacillus subtilis subsp. subtilis with (Figure 6).

Figure 6: Phylogenetic tree showing relationship of B. subtilis strain PLSSC with other Bacillus species.

The mol% G+C is an old but useful molecular taxonomy method for bacterial classification. Based on its whole genome sequence, the mol% G+C for strain PLSSC was found to be 43.58%. The average nucleotide identity (ANI) between B. subtilis strain PLSSC and B. subtilis strain KCTC 3135 was found to be ~99%. BLASTN search of the assembled PLSSC genome against the RefSeq genome database belonging to B. subtilis (taxid: 1423) showed ~99% sequence homology with the genome of B. subtilis KCTC 3135.

Comparison of B. subtilis strain PLSSC and B. subtilis KCTC 3135 genomes by BLAST Ring Image Generator (BRIG) [29] has been represented in the form of a circos plot (Figure 7). The taxonomic analyses mentioned above identify the strain as B. subtilis and reveal its closeness to B. subtilis KCTC 3135.

Figure 7: Circos plot comparison of B. subtilis strain PLSSC (size of genome is 4204670 bp) with B. subtilis strain KCTC 3135 (CP015375.1).

Antimicrobial resistance (AMR) genes

Homology search of the assembled PLSSC genome against the CARD database resulted in the identification of 717 putative antibiotic resistance genes. Identified genes were mainly involved in Defence mechanism (474), Signal transduction mechanisms; Transcription (136), General function prediction only (38), Carbohydrate transport and metabolism; Amino acid transport and metabolism; Inorganic ion transport and metabolism; General function prediction only (32), Cell wall/membrane/envelope biogenesis (10), Coenzyme transport and metabolism; Energy production and conversion (06), Replication, recombination and repair (2), Coenzyme transport and metabolism; Energy production and conversion (1), tunicamycin resistance (1), and Aminoglycoside 6-adenylyltransferase (1).

As per WHO, 2016 and EFSA, 2012 guidelines, assembled PLSSC genome was also screened for genes coding for critically important antimicrobials (CIAs) or highly important antimicrobials (HIAs). Genes for tunicamycin resistance (tmrB: DUT89_01190) and beta– lactamase (aadK: DUT89_13540) were identified in the genome. tmrB and aadK code for an ATP binding tunicamycin resistance protein and a streptomycin modifying enzyme respectively. No mobile element was identified in the flanking regions of these AMR genes which shows that these genes contribute to the intrinsic resistance and these is no risk of horizontal transfer of AMR genes.

Analysis of virulence factor genes

BLAST search against VFDB revealed that a total of 687 virulence factor proteins showed considerable homology with the assembled PLSSC genome. Further analysis with respect to the COG database suggested that these proteins were the products of nonclassical virulence factors genes and they were related to diverse COG categories. Gene mining was done to identify genes linked to Diarrheal enterotoxin bceT, Haemolytic enterotoxin operon (hbl genes – hblA, hblC, hblD), Non-haemolytic enterotoxin operon (nhe ABC genes – nheA, nheB, nheC), Cytotoxin K (cytK), Enterotoxin FM (entFM), and Emetic Toxin Cereulide (cesB). None of these genes were present in the B. subtilis strain PLSSC assembled genome.

Analysis of biogenic amine producing genes

One gene coding for amino acid decarboxylase i.e. arginine decarboxylase (DUT89_07295) was identified with 100% homology against the biogenic amine producing proteins. In the qualitative detection of activity of arginine decarboxylase gene as described by Chang et al. [30], we did not observe change in the intensity of purple colour around the colony with and without the arginine (Data not shown). Above finding confirms that the arginine decarboxylase gene is either non-functional or not expressed at a sufficient level to produce detectable amounts of biogenic amine under the tested conditions.

Analysis of genome stability

Two putative prophage regions were identified at genomic coordinates: 1176131-1209853 (33722 bp) and 2006592-2145332 (138754 bp). However, post-analysis, it was observed that both these putative prophage regions did not contain the genes required for replication/transcription/packaging (topoisomerase, replisome, DNA-binding proteins), morphogenesis (accessory – tail fiber/ whisker) or lysis (lysine). This suggests that the prophage sequences in the assembled genome were defective and non-functional.

IS finder [31] identified 12 insertion sites (IS element regions) in the assembled genome. Additionally, a search against the ACLAME database [32] revealed that 497 regions in the assembled genome showed significant hits (coverage >=50% and e-value <=1e-05). However, none of the genes coding for virulence factors had any mobile elements in their flanking regions, and therefore, do not pose any safety concerns to humans or animals. Two CRISPRs were identified from the assembled genome of B. subtilis strain PLSSC.

In silico analysis of probiotic features

It was observed that the assembled genome of B. subtilis strain PLSSC codes for important proteins which contribute to adhesion, acid and bile salt tolerance, and environmental stress resistance (Table 2).

Analysis using the antiSMASH program [33] showed that the assembled genome of B. subtilis strain PLSSC has 11 potential gene clusters responsible for the synthesis of secondary metabolites, including antimicrobial peptides, terpenes, fatty acids and others. The gene clusters are responsible for the synthesis of non-ribosomal cyclic lipopeptides: surfactin, bacillaene, fengycin, bacillibactin and bacilysin (Table S1).

Table 2: Proteins of B. subtilis strain PLSSC involved in probiotic features.