3 Biotech (2018) 8:293 https://doi.org/10.1007/s13205-018-1313-0

ORIGINAL ARTICLE

Computational identification and evolutionary analysis of toxins in Mosquitocidal Bacillus thuringiensis strain S2160-1 Panpan Liu1 · Yan Zhou2,3 · Zhongqi Wu1,3 · Hao Zhong2 · Yanjun Wei6 · Youzhi Li2 · Shenkui Liu1 · Yan Zhang6 · Xuanjun Fang1,2,3,4,5 Received: 22 March 2018 / Accepted: 18 June 2018 © Springer-Verlag GmbH Germany, part of Springer Nature 2018

Abstract Mosquitocidal Bacillus thuringiensis (Bt) strain S2160-1 was proposed to be an alternative to Bacillus thuringiensis subsp. israelensis (Bti). Discovering and validating a toxic gene by experimentation was a complex and time-consuming task, which can benefit from high-throughput sequencing analysis. In this research, we predicted and identified toxic proteins in the strain S2160-1 based on the draft whole genome sequence data. Through a local BLASP, 46 putative toxins were identified in S2160-1 genome, by searching against a customized B. thuringiensis toxin proteins database containing 653 protein or peptide sequences retrieved from public accessible resources and PCR/clone results in our laboratory (e value = 1e − 5). These putative toxins consist of 42 to 1216 amino acids. The molecular weights are ranged from 4.86 to 137.28 kDa. The isoelectric point of these candidate toxins varied from 4.3 to 10.06, and 16 out of which had a pH greater than 7.0. The analysis of tertiary structure and PFAM domain showed that 12 potential plasmid toxins may share higher similarity (9/12 QMEAN4 score > 0.3) with known Bt toxins. In addition, functional annotation indicated that these 12 potential toxins were involved in “sporulation resulting in formation of a cellular spore” and “toxin activity”. Moreover, multiple alignment and phylogenetic analysis were carried out to elucidate the evolutionary relationship among 101 known crystal or toxin proteins from public database and them with MEGA 6.0. It indicated that PS2160P2_1 and PS2160P2_153 may be potential Cry4like toxins in Bt S2160-1. This research may lay the foundation for future functional analysis of Bt S2160-1 toxin proteins to reveal their biological roles. Keywords  Bacillus thuringiensis · Crystal proteins · Phylogenetic analysis · Evolution · Genome Abbreviations CTGs Candidate toxin genes TER-CTGs Candidate toxin genes with tertiary structure

CD-CTGs Candidate toxin genes with conserved domain GO-CTGs Candidate toxin genes with GO annotations

Introduction Panpan Liu, Yan Zhou and Zhongqi Wu contributed equally to this work. Electronic supplementary material  The online version of this article (https​://doi.org/10.1007/s1320​5-018-1313-0) contains supplementary material, which is available to authorized users.

Bacillus thuringiensis is an aerobic Gram-positive sporeforming entomopathogenic bacterium. It can be isolated from a wide variety of environmental sources such as soil,

* Panpan Liu [email protected]

3



Hainan Institute of Tropical Agricultural Resources, Hainan, China

* Xuanjun Fang [email protected]

4



Institute of Life Science, Jiyang College of Zhejiang A&F University, Zhejiang, China

5



Cuixi Academy of Biotechnology, Zhejiang, China

6



School of Life Science and Technology, Harbin Institute of University, Heilongjiang, China

1



Northeast Forestry University, Harbin 150040, China



College of Life Sciences and Technology, Guangxi University, Guangxi, China

2

13

Vol.:(0123456789)

293  

Page 2 of 8

3 Biotech (2018) 8:293

water, dead insects, plant surfaces and insect feces (Federici 1999). In the unfavorable growth conditions, B. thuringiensis might generate spores, and most of the bacteria are accompanied by the production of a crystalline protein, which was called parasporal crystal or δ-endotoxin. To different insect larvae, the parasporal crystals have different pathogenicity, and some of them such as parasporins were reported to be able to kill cancer cells (Nagata et al. 2009; Palma et  al. 2014; Poornima et  al. 2010; Schnepf et  al. 1998). In the families of the delta endotoxins produced by B. thuringiensis, Cry toxins and Cyt toxins were the most widely studied families (Jouzani et al. 2017). In past few decades, more than 700 genes encoding Cry crystal proteins were identified and characterized in Bt Toxin Nomenclature (http://www.lifesc​ i.sussex​ .ac.uk/Home/Neil_Crickm ​ ore/Bt/) (Crickmore and Feitelson 2016). As the insecticidal crystal protein had specific insecticidal activity to target insects but safe to non-target organisms, it was widely used for controlling insect pests in agriculture and forestry as a biological control agent. However, long term use of the same Bt insecticide or single cropping of insect-resistant transgenic crop will lead to insect resistance to Bt toxin protein, therefore, it is very urgent to constantly explore new Bt strains and novel toxins for the pest control (Bellile and Vonesh 2016; Tabashnik et al. 2013). Cry1Ab, which was identified and cloned by molecular biological approaches from B. thuringiensis strain SY49-1, was reported with insecticidal activity, and the insecticidal activity had been improved by the technique of site-directed mutagenesis. Moreover, a new B. thuringiensis with multiple insecticidal effects showed previously was constructed by means of genetic engineering, binding transfer, etc. (Azizoglu et al. 2016; Zhang et al. 2012). However, these methods of driftlessly mining and modifying new toxin molecules need doing a lot of experiments and consuming a large number of manpower and material Fig. 1  The workflow of computational identification and evolution analysis for toxic proteins in Bt S2160-1

resources, which is time-consuming, costly and inefficient. In an improvement over conventional experiment screening, next generation sequencing (NGS) has been used for identification of new toxins from Bt strains. To narrow the scope of candidate toxins by screening test and improve the identification efficiency of toxins, here, a bioinformatics method was proposed to computationally identify toxins in whole genome based on sequence similarity and conserved functional domains of toxins. Bt S2160-1 with acute mosquitocidal activity was isolated from soil samples from Dawangling Forest Nature Reserve (Guangxi, China) and proposed as an alternative to Bti for mosquitocidal products (Zhang et al. 2012). In 2013, the whole genome sequencing of Bt S2160-1 have been done and its draft genome was assembled that contains three replicas, a nucleoid genome in length of 5.46 Mb and two plasmid genomes with the lengths of 544 and 309 kb (Zhong et al. 2013). Cry4Cb3 with mosquitocidal activity against Culex quinquefasciatus was identified with mass spectrometry (Zhang et al. 2014). Nevertheless, there have been no systematic studies of the toxic protein of the Bt strain S21601. With the completion of whole genome sequencing of Bt S2160-1, we had opportunity to investigate a genome-wide analysis of crystal proteins. In this study, an in silico search of Bt 2160-1 genome database was conducted to identify crystal proteins. We applied multiple tools and programs to conduct in-depth analysis of each toxin proteins, including the physiological and biochemical properties, sequence structures, functional domains. Then, the GO function annotation was proceeded to study the potential functions of crystal toxins. Furthermore, a thorough comparative analysis of Bt 2160-1 toxin proteins to those from other Bt strains was performed (Fig. 1). Our results presented here may provide a subset of potential candidate toxins for

Completed work Bt Strain S2160-1 whole genome sequencing Bt Strain S2160-1 draft genome Predicted peptide sequences

Candidate Toxin Genes (CTGs)

Physical and chemical properties analysis Phylogenetic Analysis Protein tertiary structure (TER-CTGs) prediction analysis

Conserved function domains (CD-CTGs) identification

Function annotation with GO

13

Screened Toxin Genes (STGs)

known Bt toxins from UniprotKB

Page 3 of 8  293

3 Biotech (2018) 8:293

future engineering modifications of mosquitocidal activity characteristics in Bt 2160-1.

Materials and methods Sequence retrieval and construction of the local protein database for toxins Bacillus thuringiensis Toxin Nomenclature (http://www. lifes​ci.susse​x.ac.uk/home/Neil_Crick​more/Bt/intro​.html) (Crickmore and Feitelson 2016): The Bacillus thuringiensis delta endotoxin nomenclature committee was set up in 1993 to update the nomenclature originally devised in 1989 by Hofte and Whitele. The corresponding protein sequence of delta endotoxins within nomenclature was downloaded from NCBI according to the GI and accession numbers (up to July, 2016). In addition, sequences cloned by our lab that have not submitted to public database were also included. Finally, 653 toxin proteins were selected to create a local toxin protein database. To preliminarily filter out the potential toxins, the peptide sequences of Bt S2160-1 were used as query sequences to search against the local database via local BLASTP (e value = 1e − 5).

Sequence analysis of potential toxins in Bt strain S2160‑1 To make a preliminary judgment on whether a candidate protein meets the physicochemical parameters of toxic proteins, the theoretical isoelectric point (pI), molecular weight (MW) and amino acid composition of these candidate toxins in Bt strain S2160-1 were calculated by ProtParam (http://web. expasy​ .org/protpa​ ram/) (Wilkins et al. 1999). SOPMA (https​ ://npsa-prabi​.ibcp.fr/cgi-bin/npsa_autom​at.pl?page=npsa_ sopma.​ html) was used to predict the secondary protein structure of these candidate toxins in S2160-1 with default parameters: Window width: 17, Similarity threshold: 8, Number of states: 4. Then, protein sequences of these candidate toxins were uploaded to SWISS-MODEL (http://swissm ​ odel.​ expas​ y.org/) (Biasini et al. 2014; Kiefer et al. 2009) to build the protein tertiary structure homology modelling. To explore the domain organization of these candidate toxins, PFAM (http://pfam.xfam.org/) (Finn et al. 2016) was performed. Moreover, the functional annotation of these potential toxins was performed using Blast2GO program (Conesa and Gotz 2008), including three parts: biological process, molecular function and cellular component.

Phylogenetic and sequence alignment analysis of hypothetical toxins in Bt S2160‑1 Multiple sequence alignments of hypothetical toxin protein sequences of Bt S2160-1 and other known Bt toxins

from UniprotKB were performed using ClustalW tool of MEGA 6.0 (Molecular Evolutionary Genetics Analysis) with default parameters. Phylogenetic tree was constructed using maximum likelihood method with 100 replicates of bootstrap analysis in MEGA 6.0 (Tamura et al. 2013). The phylogenetic tree visualization and later landscaping is made with iTOL (Interactive Tree of Life) online tool (Letunic and Bork 2016).

Results Computational identification of Bt strain S2160‑1 candidate toxins A customized local protein database was constructed with 621 sequences manually maintained by Crickmore et al., and 32 genes cloned by our research team but not registered in public databases (up to 2016-07). BLASTP (version 2.2.28) was performed locally to identify potential toxic proteins (e value = 0.00001) in Bt strain S2160-1. Each blast hit was assigned to the corresponding group incorporating four ranks as described by Crickmore et al., in which the boundaries represent approximately 95, 78, and 45% sequence identity (Crickmore et al. 1998). Finally, we identified 46 candidate toxin genes (marked as CTGs), consisting of 6 potential vip, 37 cry and 3 cry/cyt genes (Supplementary Materials Table S1). For the sake of convenience, each gene was named as following: If genes located on circular chromosome, they were named like nS2160_gene_id_2 and the corresponding encoding protein was named like NS2160_2. Otherwise, genes mapped on plasmids were named as pS2160p1_gene_id_2 and the corresponding proteins were named as PS2160P1_2.

Computation of physical and chemical parameters for the candidate toxic proteins As described previously, we collected 653 toxins, which were currently known or kept only in our lab. Results of ProtParam indicated that the conserved Bt toxin proteins held the following physical and chemical properties: (1) count of amino acids ranged from 78 to 1522, (2) MW ranged from 8.59 to 170.47 kDa, and (3) pI ranged from 4.07 to 9.79. To determine whether these 46 candidate toxic proteins had the similar features as the conserved Bt toxin proteins, we computed the physical and chemical parameters for all of them. It showed that protein sequences of these CTGs contained 42 (NS2160_814) ~ 1216 (PS2160P2_1) amino acids (aa) (Supplementary Materials Table S2). The MW of them ranged from 4.86 to 137.28 kDa, and the pI value ranged from 4.3 (PS2160P1_40) to 10.06 (PS2160P2_241).

13

293  

Page 4 of 8

It seemed that the physical and chemical properties of these CTGs were similar to that of the conserved Bt toxin proteins.

Structure and conserved motif analysis of the candidate toxic proteins To investigate whether proteins coded by these 46 CTGs in Bt S2160-1 had homological model with known Bt insecticidal crystal protein, SWISS-MODEL (https​://www.swiss​ model​.expas​y.org/) (Biasini et al. 2014; Kiefer et al. 2009) was used to predict the protein tertiary structures. It showed that 14 of 46 candidate proteins had at least one matching result with the three-tier models of known Bt insecticidal crystal proteins (Supplementary Materials Table S3). The genes coding these 14 proteins were designed as candidate toxin genes with tertiary structure (marked as TER-CTGs). Typically, the active conformation of Bt δ-endotoxins was composed of three distinct structural domains: Endotoxin_N, Endotoxin_M and Delta_endotoxin_C (Bravo et al. 2007; Galitsky et al. 2001; Palma et al. 2014), which play important roles in entomopathogenicity, membrane insertion and pore formation, receptor binding, etc. Some crystal proteins do not have typical endotoxin domains, but contain RICIN or Toxin_10 domain (Itsko et al. 2005; Kelker et al. 2014). The distributions of typical endotoxin domains in these 46 candidate toxin proteins were investigated with PFAM (http://pfam.xfam.org/) (Finn et al. 2016). Fortunately, we found that (1) NS2160_5341, PS2160P1_21, PS2160P1_40, PS2160P1_49 and PS2160P2_146 had δ_endotoxin_C domain, (2) PS2160P1_2, PS2160P1_20, PS2160P1_39, PS2160P1_50, PS2160P1_54, PS2160P2_1, PS2160P2_134, PS2160P2_147, PS2160P2_157, and PS2160P2_246 had Endotoxin_N, Endotoxin_M and Delta_ endotoxin_C domains (Supplementary Materials Table S3). As far as is known, most toxic genes if not all Cry and Cyt protein-coding genes were mapped on plasmids present in Bt, not on the chromosome. In our results, all of them except NS2160_5341 were located on plasmids. So, we marked them as candidate toxin genes with conserved domains (CD-CTGs).

GO functional annotation for CTGs The sequence diversity of Bt strain S2160-1 toxins indicated that these toxin proteins may be involved in multiple biological processes related to pathogen. To understand the biological processes associated with Bt strain S2160-1 toxins, we performed GO annotation of the candidate toxins. Finally, 31 of 46 CTGs were annotated to GO terms (marked as GO_ CTGs). The results are shown in Supplementary Materials Table S4. The majority of toxic protein-coding CTGs located on plasmids were involved in binding (25 CTGs), cellular process (19 CTGs), multi-organism process (17 CTGs),

13

3 Biotech (2018) 8:293

interspecies interaction between organisms (17 CTGs), sporulation (14 CTGs), development process (14 CTGs) and signaling receptor binding (12 CTGs). Meanwhile, nine plasmid proteins and four chromosome proteins participated in metabolic process, organic substance metabolic process, nitrogen substance metabolic process, cellular metabolic process and primary metabolic process. From which, it was reasonable to speculate that these 13 proteins may not true toxins in Bt strain S2160-1. In addition, three proteins (NS2160_3788, NS2160_3497, NS2160-3703) located on S2160-1 chromosome functioned in biosynthetic process, small molecule metabolic process and oxidation–reduction process. Based on the GO functional annotation analysis, we found that CTGs of Bt strain S2160-1 could function in the development process of sporulation, signaling receptor binding and interspecies interaction between organisms. Summarily, we identified 46 potential toxins in Bt strain S2160-1 via BLASTP. Combining the analysis of tertiary structure, conserved domains and GO annotation of these candidate toxins, 12 plasmid genes of 46 CTGs from S2160-1 were more likely to be a Bt toxin (Fig. 2a). From the functional analysis of these CTGs, the intersection of CDS_CTGs, TER_CTGs and GO_CTGs was assigned into GO:0090729 (MF: toxin activity), GO:005102 (MF: signaling receptor binding), GO:0009405 (BP: pathogenesis) and GO:030435 (BP: sporulation resulting in formation of a cellular spore) (Fig. 2b). Moreover, PS2160P1_49 and PS2160P1_50 were annotated as mosquitocidal toxin-like protein, which was consistent with mosquitocidal activity of Bt S2160-1 (Fang and Zhang 2012) (Table 1).

Phylogenetic analysis of the toxin genes in 12 Bt strains To investigate the phylogenetic relationships of toxins between Bt S2160-1 and that of other Bt strains, a phylogenetic tree was constructed with maximum likelihood (ML) method for 100 bootstraps in MEGA 6.0 (Tamura et al. 2013). The phylogenetic tree included the above 12 computational predicted toxic proteins from Bt S2160-1 and 101 known toxins from UniprotKB (Fig.  3). From the evolution tree, most of the computational predicted toxic proteins (9/12) from S2160-1 had higher similarity (bootstrap value > 0.68) with toxins from Bacillus thuringiensis subsp. Israelensis (BACTI), which was known as mosquitocidal model Bt strain that could produce specific insecticidal toxins during sporulation and had been widely used in the larval mosquito controls (Guan et al. 2012, 2014). On the basis of phylogenetic analysis, we found that PS2160P2_1, PS2160P2_153 and PS2160P2_147 had similarity with Cry4 class crystal proteins. In detail, PS2160P2_1 showed high identities of less than 45% to the best hit Cry4Ba_BACTI, PS2160P2_147 had high

Page 5 of 8  293

3 Biotech (2018) 8:293 Fig. 2  The Venn diagrams of CTGs, TER-CTGs, CD-CGTs and GO-CTGs and the network of genes and GO terms. a CTGs represent the 46 candidate toxin genes identified with BLAST. TER-CTGs indicate those 14 toxin gene-coding proteins with tertiary structure. CD-CTGs represented the 14 toxin genecoding proteins with conserved endotoxin domains. GO-CTGs indicated 31 genes with GO annotation. b The network of genes and GO terms. The circle nodes in the network represent CTGs, and the square nodes represent GO terms. The color of square nodes was colored with the corresponding category

a

CDS_CTGs

TER_CTGS CTGs

0

0

0

2

GO_CTGs

GO_CTGs

15

b

TER_CTGs

1

CTGs

Biological processes related to pathogen or toxicity

GO:0006508 interacts with PS2160P2_146 PS2160P2_1 NS2160_2305 interacts with PS2160P1_40 interacts with with GO:0004222 PS2160P2_147 interacts interacts with interacts with PS2160P1_2 interacts with with PS2160P2_155 interacts withinteracts interacts with GO:0004803 interacts with interacts with interacts with interacts with interacts with interacts with interacts with interacts with interacts with PS2160P2_247 GO:0005102 interacts with interacts with PS2160P2_153 interactsinteracts with interactswith with interacts with PS2160P1_210 interacts with interacts with interacts with interacts with with interacts with interacts GO:0016021 interacts with interacts with interacts with GO:0090729 interacts with interacts with interacts with PS2160P1_119 interacts with interactsinteracts with PS2160P1_49 interacts with interacts with with PS2160P2_134 interacts interacts with interacts with with GO:0006313 interacts with interacts interactswith with interacts interacts with interacts with interacts with with GO:0030435 interacts with interacts PS2160P1_33 interactswith with interacts with with interacts with interacts interacts with with PS2160P1_20 interacts interacts with GO:0009405 with withwith GO:0003677 interactsinteracts with interacts interacts PS2160P1_50 withinteracts interactsinteracts with interacts interacts with with with GO:0016020 interacts with with with PS2160P2_208 interacts interacts interacts with PS2160P1_54 interacts with PS2160P1_39 NS2160_5338 PS2160P1_21 NS2160_5341 GO:0043565 GO:0003855 GO:0003676 GO:0008081 GO:0055114 interacts with with interacts with with interacts GO:0008270 interacts with PS2160P2_148 interacts GO:0016829 PS2160P1_24 interacts with GO:0009073 NS2160_3497 interacts with interacts with interacts with interacts with interacts with PS2160P1_8 NS2160_3703 GO:0008835interacts GO:0009231 GO:0015074 with interacts with interacts with interacts interacts with with interacts with interacts with interacts with NS2160_3788 GO:0009423 interacts with with interacts PS2160P2_135 PS2160P2_171 PS2160P1_30 GO:0050661 GO:0006629 GO:0008703 Molecular Function Intersection of CDS TER and GO_CTGs

Table 1  The general features of 12 toxin genes in Bt strain S2160-1

CDS_CTGs

13

Cellular component Intersection of CDS and GO_CTGs

Biological Process Intersection of TER and GO_CTGs

GO_CTGs

Protein name

Number of aa

Genome location

Description with GO

PS2160P1_2 PS2160P1_20 PS2160P1_21 PS2160P1_39 PS2160P1_40 PS2160P1_49 PS2160P1_50 PS2160P1_54 PS2160P2_1 PS2160P2_134 PS2160P2_147 PS2160P2_153

699 525 69 601 96 117 575 688 1216 802 552 1157

p01:2258-4357 p01:23917-25491 p01:25575-25784 p01:42105-43907 p01:43985-44275 p01:50832-51185 p01:51263-52987 p01:55475-57541 p02:216-3866 p02:114830-117238 p02:129794-131449 p02:136654-140127

Cry-like protein Cry30-like protein Cry30-like protein Pesticidal crystal protein Cry53Ab1-like protein Mosquitocidal toxin-like protein Mosquitocidal toxin-like protein Pesticidal crystal protein Cry-like protein Cry-like protein Cry50-like protein Pesticidal crystal protein

13

3 Biotech (2018) 8:293

known toxins

AC

CU

H

CU

F

CT

U

CT

Cry4Aa BACTI

0.52 1

Spo0A B ACTU SP2G BACT K Cry1 5Aa 0.0 BAC 1 UT 0.96 PS2 160P 1 40 PS2 1 60P 0.9 1 49 2 RP SX 0 .7 BA R C 4 P2 TK 0.9 8B 9 AC RP TK 35 B Cr AC y1 TK 3A Cr aB y1 AC 4 Cr Aa TU y2 BA 1A Cr CT y5 a S B Ba AC BA TU CT U

BA

BA

0.03

BA

24

1

Aa B

0.33

0.56 1

9

7 0. 55

1

85

0. 25

1

TM

AC

CT

aB

BA

50

T1

aB

AC

aB

TF

AC

BA y3 CT Ba Cr U BA y3 Ca Cr CT y3 O B AC Aa TK BA CT T

TF

G

CT

y1

Cb

Bb

BA C

CU

BA

Cr y1

8A

P1

BA

Cr

Aa

y8

TE Ja G BA Cry CT 1Jb U BA Cry CT 1Ba U BA CT Cry E 1Ba BAC TK Cry1 Be B 1 ACT U Cry1 0.99 Bc B ACTM Cry1B b BAC 1 TU 1 Cry1Bd B ACTZ 1 Cry1Ka BACT M

CT A CT A

y2

y3

TI

160

y8

Cr

BA

Cr

Cr

Cr

BAC

TU

39

Ba

Cr

y1 C

d BAC

60P1

y8

Cr

U

TP

TO

TA

AC

AC

aB

aB

9E

Aa

y9

Cry

Cr

9C

Cry

ACT

BAC

Aa B

9Da

Cry7

Cry

UA

K ACU

Ab B

Cry7

ACTU

b BAC

Cry1Id B

Cry1Ib BACTE

Cry1Ia BACTK

Cry7A

Cr

9

26A

U

BACTK

TD C M BA a CT 3A BA TP ry Aa C AC B K y3 Cr CU Ca BA K

8

1

3 0.9 0.55

BA

aS

Cry

0.9

0.

0.9

C

Cd

1

0.75

Ca

PS2

1

BA

87

2

1F b

Hpr

0.9

Fa

Cry2A PS21

1

Cr y1

Cry2Ab

1

0.

0.

5

0.96

ry

2

5 0.6

0.4

Cry2Aa BACTX Cry2Ac BACT

0.6

0.85

1

C

65

0.97

y1

Cry2Aa BACTK

0.45 0.97

0.

0.62

ACTK 0.81 Cry1Ab B 1 CTB A B b Cry1A CTU Af BA Cry1 1 CTK A Ac B A Cry1 ACT B 1Eb X Cry CT BA 1 Ea 1 M 1 y CT Cr 0.9 BA b U T 1H C y r 1 C BA U Ha CT y1 BA Cr A b CT D y1 U BA Cr CT Da 1 A Z y B Cr a CT G A B y1 b Cr G y1 Cr

Cr

V

Cry11Bb BACT

0.98

0.98 0.93

0.45

1

A Cry1Ab BACT

U

CT

BA

1

1

Cry1Ae BACTL

2

0.66

Cry1Ic BACTU

Cry1Ad BACTA

1

Cry1Aa BA CTK Cry1Ag BACTU

6

0.89 0.85 1 0.85 0.27

0.99

TE

0.5

0.4

0.3

a BAC

1

Cry1A

1

ACT S ACTA

Aa B

1

Cry1

TU C

BA

C

0.

U

TK

Aa

D CU Ac BA D y5 b Cr CU A y5 BA TY a Cr A AC y5 1B Cr a TI t2A AC 97 65 1 .82 Cy aB 0. 0. J B 0 t2 CT Cy BA Bb W 2 t CT Cy BA 9 9 1Ba TI 0.8 0.9 Cyt C BA 1Aa 1 M Cyt C 0.1 A T 0.9 Aa B Cyt1 V 8 T 9 C . A 0 .78 bB 0 .45 Cyt1A 0 CTI A B a Cry11A 0.59 BACTJ Cry11Ba 1 ry

0. 89 0. 6

NB

S2160-1 toxins

a

BA

0.65

Aa

a

0A

Ba

y6

7A

0.98

Cr

To x

24 ry

y6

25

y2

Cr

Cry4 class toxins

Cry1

y2

Cr

Cry2 class toxins

C

Cr

Cyt class toxins

1

y Cr

Cry11 class toxins

PS2160P2 153

Cry class toxins

Cry4Ba BACTI

147 PS2160P2 0P2 1 PS216 54 60P1 0.88 86 PS21 0. 1 21 0.63 160 P PS2 20 P1 9 160 0 .9 TI PS2 AC aB 10A 12 9 01 Cry 0P 216 0.9 0. UH PS AC J aB T C 9B A y1 4 aB Cr 03 13 9A 0. y1 P2 TU Cr 60 21 AC J B T PS C a 2A BA TJ y2 C Aa Cr BA Aa

Colored ranges

Cry4AA BACT K

Page 6 of 8

1

293  

Fig. 3  The maximum likelihood phylogenetic tree of toxin genes in Bt strains. 12 CTGs of Bt S2160-1 and 101 known toxins from UniProtKB database were used to construct the tree with bootstrap

repeated 100 times, and the bootstrap value larger than 0.50 was shown in the figure. Proteins from S2160-1 were marked with red circle. Each known toxin was marked with gray circle

identities of 53% to its best hit Cry4Ba_BACTI and PS2160-1_153 had high identities of 48% with its best hit Cry4Aa_BACTI. Meanwhile, PS2160P1_2, PS2160P1_20, PS2160P1_21, PS2160P1_54 had closed evolutionary relationship with Cry10 (bootstrap value > 50). In summary, the construction of phylogenetic tree could better reflect the relationship among the toxin genes in Bt strains that could be helpful to classify and name the toxin genes as well as to infer the function of toxins in Bt S2160-1 strain.

Discussion

13

B. thuringiensis (Bt) is a natural typical insecticidal bacterium. Using Bt bacteria as biological pesticides or developing transgenic Bt crops is becoming the main directions of research and development in agrochemical and seed industry. Bt always performed substantial variations in insecticidal toxicity. Therefore, how to quickly identify new Bt strains and how to discover new Bt toxic proteins are the basis for establishing competitive advantages in these industries.

Page 7 of 8  293

3 Biotech (2018) 8:293

The development of genome sequencing technology, especially the next generation sequencing, provides new means for discovering Bt toxins in large scale. In this study, we employed the data of next generation sequencing (NGS) data and bioinformatics approaches to computationally discover and identify Bt toxin proteins in Bt strain S2160-1, which demonstrated that it was an efficient strategy comparing with the previous traditional PCR method. In theory, NGS could produce a large number of sequences quickly and indiscriminately. On the other hand, the standard computational process could identify all types of B. thuringiensis toxin proteins, including Cry proteins, Cyt proteins and Vip proteins. Thus, it was expected that computational method was more efficient and time-saving in identifying Bt toxins than traditional PCR-based systems. In fact, we identified 46 candidate toxins for Bt strain S2160-1. The tertiary structure analysis with SWISS-MODEL showed that 15 of 46 proteins (TER-CTGs) encoded by the CTGs had a matching result with the three-tier models of known Bt insecticidal crystal protein. Meanwhile, 12 proteins (CDCTGs) among these 15 TER-CTG-encoded proteins were located on circular plasmids, which contained at least one of delta-endotoxin domains. The functional annotation showed that they might be involved in the biological processes or molecular function related to Bt pathogen (Fig. 2b). The good news was that four toxins of the insecticide genes Cry30Ea1 (NCBI Acc No. ACC95445), Cry30Ga2 (NCBI Acc No. HQ638217), Cry50Ba1 (NCBI Acc No. GU446675) and Cry54Ba1 (NCBI Acc No. GU446677) previously cloned in Bt S2160-1 were confirmed by our results (Zhang et  al. 2012) (PS2160P1_54 corresponded to Cry30Ea1, PS2160P1_20 corresponded to Cry30Ga2, PS2160P2_147 corresponded to Cry50Ba1, and PS2160P1_2 corresponded to Cry54Ba1). Additionally, PCR approach had identified three distinguishing bands of parasporin proteins with 140, 130 and 30 kDa in size but remained unknown, which was perfectly matched with this study that PS2160P2_1 was predicted to be with 137.28 kDa and PS2160P2_153 was predicted to be with 130.57 kDa. Combining the results of BLASTP and phylogenetic tree, it was reasonable to infer that PS2160P2_1 and PS2160P2_153 may be potential Cry4-like toxins in Bt S2160-1. It was obvious that our research might facilitate to detect those Cry toxins with three-domain through bioinformatics methods. Remarkably, our research also identified those candidate toxins without typical three domains (Kirouac et al. 2002; Vachon et al. 2012). For example, the protein coded by pS2160p1_gene_id_8 did not contain typical delta endotoxin domains, whereas contained RICIN domain as similar as Cry35 class crystal proteins and Cyt1Ca, and with an active tertiary structure prediction result. Previous researches reported that the Bt toxins Cyt1Ca contained RICIN domain, and had 50% similarity to the β-trefoil modules which was

found in various toxins containing RICIN domains, such as the mosquitocidal toxin protein Mtx1 from Bacillus sphaericus (Berry et al. 2002; Manasherob et al. 2006). Therefore, it was reasonable to speculate that the protein coded by pS2160p1_gene_id_8 may be a member of two-domain fusion toxin as Cyt1Ca. In conclusion, this computational research based on NGS was proposed to mine cry genes from B. thuringiensis. The method was able to mine toxins in an efficient and fast way. As a result, 46 candidate toxins were identified from S21601. The GO function annotation revealed that 23 of them were annotated to known pesticidal crystal proteins. Among them, 12 toxins with typical three domains can be identified accurately. In the end, 12 proteins with all features were chosen to investigate their evolutionary relationship with known toxins. The above results proved the efficiency of this method in mining toxin genes. In fact, the mining of novel sequences is associated with the selection of known reference sequence. It was important to note that the quality of sequence and assembled genome affects the final prediction results, so a choice would need to be made between the cost and efficiency of sequencing. And this method could mine all classes of toxin proteins based on sequence similarity. This research might make key contributes to further study on the generic features of B. thuringiensis S2160-1 toxins. Acknowledgements  This work was part of the project “China National Bt Strains Resource Initiative (BtSRI)” funded by Hainan Institute of Tropical Agricultural Resources in Hainan province (HITAR). All the intellectual property rights are owned by HITAR in Hainan province.

Compliance with ethical standards  Conflict of interest  No conflict of interest declared.

References Azizoglu U, Ayvaz A, Yilmaz S, Karaborklu S, Temizgul R (2016) Expression of cry1Ab gene from a novel Bacillus thuringiensis strain SY49-1 active on pest insects. Braz J Microbiol [Publ Braz Soc Microbiol] 47:597–602. https​: //doi.org/10.1016/j. bjm.2016.04.011 Bellile KG, Vonesh JR (2016) Bioinsecticide and leaf litter combination increases oviposition and reduces adult recruitment to create an effective ovitrap for Culex mosquitoes. J Vector Ecol J Soc Vector Ecol 41:123–127. https​://doi.org/10.1111/jvec.12203​ Berry C et al (2002) Complete sequence and organization of pBtoxis, the toxin-coding plasmid of Bacillus thuringiensis subsp. israelensis. Appl Environ Microbiol 68:5082–5095 Biasini M et al (2014) SWISS-MODEL: modelling protein tertiary and quaternary structure using evolutionary information. Nucl Acids Res 42:W252–W258. https​://doi.org/10.1093/nar/gku34​0 Bravo A, Gill SS, Soberon M (2007) Mode of action of Bacillus thuringiensis cry and Cyt toxins and their potential for insect control. Toxicon Off J Int Soc Toxinol 49:423–435. https​://doi. org/10.1016/j.toxic​on.2006.11.022

13

293  

Page 8 of 8

Conesa A, Gotz S (2008) Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int J Plant Genom 2008:619832. https​://doi.org/10.1155/2008/61983​2 Crickmore NZD, Feitelson J (2016) Bacillus thuringiensis toxin nomenclature http://www.lifes​ci.susse​x.ac.uk/Home/Neil_Crick​ more/Bt/ Crickmore N et al (1998) Revision of the nomenclature for the Bacillus thuringiensis pesticidal crystal proteins. Microbiol Mol Biol Rev 62:807–813 Fang XJ, Zhang Z (2012) Bt S2160-1, a high mosquitocidal strain of Bacillus thuringiensis. Bt Res 3:29–32 doi. https ​ : //doi. org/10.5376/bt.2012.03.0005 Federici BA (1999) Bacillus thuringiensis. In: Bellows TS, Gordh G, Fisher TW (eds) Handbook of biological control. Academic Press, San Diego: 517–548 Finn RD et al (2016) The Pfam protein families database: towards a more sustainable future. Nucl Acids Res 44:D279-285 https:​ //doi. org/10.1093/nar/gkv13​44 Galitsky N, Cody V, Wojtczak A, Ghosh D, Luft JR, Pangborn W, English L (2001) Structure of the insecticidal bacterial delta-endotoxin Cry3Bb1 of Bacillus thuringiensis Acta crystallographica Section D. Biol Crystallogr 57:1101–1109 Guan P et al (2012) Complete genome sequence of Bacillus thuringiensis serovar Sichuansis strain MC28. J Bacteriol 194:6975. https​:// doi.org/10.1128/JB.01861​-12 Guan P et al (2014) Bacillus thuringiensis subsp. sichuansis strain MC28 produces a novel crystal protein with activity against Culex quinquefasciatus larvae. World J Microbiol Biotechnol 30:1417– 1421. https​://doi.org/10.1007/s1127​4-013-1548-1 Itsko M, Manasherob R, Zaritsky A (2005) Partial restoration of antibacterial activity of the protein encoded by a cryptic open reading frame (cyt1Ca) from Bacillus thuringiensis subsp. israelensis by site-directed mutagenesis. J Bacteriol 187:6379–6385. https:​ //doi. org/10.1128/JB.187.18.6379-6385.2005 Jouzani GS, Valijanian E, Sharafi R (2017) Bacillus thuringiensis: a successful insecticide with new environmental features and tidings. Appl Microbiol Biotechnol 101:2691–2711. https​://doi. org/10.1007/s0025​3-017-8175-y Kelker MS et al (2014) Structural and biophysical characterization of Bacillus thuringiensis insecticidal proteins Cry34Ab1 and Cry35Ab. PLoS One 9(1):e112555. https​://doi.org/10.1371/journ​ al.pone.01125​55 Kiefer F, Arnold K, Kunzli M, Bordoli L, Schwede T (2009) The SWISS-MODEL repository and associated resources. Nucl acids Res 37:D387–D392. https​://doi.org/10.1093/nar/gkn75​0 Kirouac M, Vachon V, Noel JF, Girard F, Schwartz JL, Laprade R (2002) Amino acid and divalent ion permeability of the pores

13

3 Biotech (2018) 8:293 formed by the Bacillus thuringiensis toxins Cry1Aa and Cry1Ac in insect midgut brush border membrane vesicles. Biochimica et Biophysica Acta 1561:171–179 Letunic I, Bork P (2016) Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees Nucl Acids Res 44:W242-245 https:​ //doi.org/10.1093/nar/gkw290​ Manasherob R, Itsko M, Sela-Baranes N, Ben-Dov E, Berry C, Cohen S, Zaritsky A (2006) Cyt1Ca from Bacillus thuringiensis subsp. israelensis: production in Escherichia coli and comparison of its biological activities with those of other Cyt-like proteins. Microbiology 152:2651–2659. https​://doi.org/10.1099/mic.0.28981​-0 Nagata H et al (2009) Identification of the binding domain of Streptococcus oralis glyceraldehyde-3-phosphate dehydrogenase for Porphyromonas gingivalis major fimbriae. Infect Immun 77:5130– 5138. https​://doi.org/10.1128/IAI.00439​-09 Palma L, Munoz D, Berry C, Murillo J, Caballero P (2014) Bacillus thuringiensis toxins: an overview of their biocidal activity. Toxins 6:3296–3325. https​://doi.org/10.3390/toxin​s6123​296 Poornima K, Selvanayagam P, Shenbagarathai R (2010) Identification of native Bacillus thuringiensis strain from South India having specific cytocidal activity against cancer cells. J Appl Microbiol 109:348–354. https​://doi.org/10.1111/j.1365-2672.2010.04697​.x Schnepf E et al (1998) Bacillus thuringiensis and its pesticidal crystal proteins. Microbiol Mol Biol Rev MMBR 62:775–806 Tabashnik BE, Brevault T, Carriere Y (2013) Insect resistance to Bt crops: lessons from the first billion acres. Nat Biotechnol 31:510– 521. https​://doi.org/10.1038/nbt.2597 Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol 30:2725–2729. https​://doi.org/10.1093/molbe​v/mst19​7 Vachon V, Laprade R, Schwartz JL (2012) Current models of the mode of action of Bacillus thuringiensis insecticidal crystal proteins: a critical review. J Invertebr Pathol 111:1–12. https​://doi. org/10.1016/j.jip.2012.05.001 Wilkins MR, Gasteiger E, Bairoch A, Sanchez JC, Williams KL, Appel RD, Hochstrasser DF (1999) Protein identification and analysis tools in the ExPASy. server. Methods Mol Biol 112:531–552 Zhang W et al (2012) Characterization of a new highly mosquitocidal isolate of Bacillus thuringiensis—an alternative to Bti? J Invertebr Pathol 109:217–222. https​://doi.org/10.1016/j.jip.2011.11.003 Zhang W et al (2014) Identification of a mosquitocidal toxin from Bacillus thuringiensis using mass spectrometry. World J Microbiol Biotechnol 30:3273–3277. https​://doi.org/10.1007/s1127​ 4-014-1744-7 Zhong HWY, Zhou Y, Huang FY, Zhang WF, Zhang Y (2013) Draft genome sequence of Bacillus thuringiensis strain S2160-1 with high mosquitocidal activity. Bt Res 4:29–31