residues implicated in the formation K⁺ and Na⁺ ion channels effector toxins

Sequences preparation and highly similar sequences elimination

Key words "K⁺, Na⁺, channel, scorpion toxin" were used to search the NCBI database [25] which is linked to swissprot, pdb and embl among other databases, for existing K⁺ and Na⁺ channel effectors (toxin sequences). Whole length and fragment sequences were included to insure maximum coverage of these toxins with sequence information.

All sequences were gathered according to their activities. Sequences that belong to the same group were put in the same file in FASTA format. Three groups were obtained: the first contains toxins active on the Na⁺ channels; the second contains toxins active on the K⁺ channels, and the third contains birtoxin-like peptides. A perl script was used to conduct batch PSI-BLAST [26] against the nr (non-redundant protein) database to search for similarity and to characterize these toxins.

Multiple sequence alignment

Multiple sequences alignment was performed using ClustalW [27] version 2 [28]. Multiple alignments of sequences for each group of toxins and for all sequences pooled together, were carried out in order to characterize potential conserved and variable areas for each sequence partition as well as consensus sequences (centroides). To date, ClustalW is still the most popular alignment tool. However, recent methods, in some cases, offer significantly better alignment quality. Thus, this same stage was carried out by another multiple alignments software: CHAOS & DIALIGN [29], which is based on a combination of local and global alignment methodologies. The obtained results were compared with the former ones.

In order to keep random noise at its lowest level, all qualified sequences contained three or four S-S pattern. Based on these results, highly similar sequences were removed to reduce potential bias on the motifs search (motifs conserved in each group). A subset of the originally found sequences, judged representative of the K⁺ and/or Na⁺ ion channel effectors, were used as training dataset. This subset reflects a trade-off between sequences that cover most of scorpion species, and sequences that are centroide in the multiple alignments. However, for the birtoxin-like group, all sequences were sampled.

Motifs search

Using the training dataset composed of the selected toxins representative of the ion channels blockers for K⁺ and Na⁺, GIBBS [30, 31] version 2.05 and MEME [32, 33] version 3.5.4 which is part of the Meta-MEME [34] package, were used to dig out conserved motif information, for all considered sequence partitions.

We used GIBBS to discover the conserved motif(s) (with pre-fixed length, which can contain(s) gaps) specific to each group of ion channel effectors. This same stage was performed by MEME, which determines un-gapped multilevel motifs (with no pre-fixed length). The most significant results were compared and retained.

Database search for new K⁺ and Na⁺ channel effector toxins

The output of these programs (GIBBS and MEME) was thoroughly investigated, and then used as input in the next step of the analysis.

A profile using these motifs was systematically formed, and used to search for sequences with this profile. For that purpose, MAST, which is also part of the Meta-MEME package, was used to search for new K⁺ and Na⁺ channel effector sequences (test dataset) against nr database. The e-value for a qualified sequence was set to 1. A perl script was used to eliminate replications (identical to the sequences of the training dataset).

Interproscan

InterPro [35] is an integrated resource for protein families, domains and functional sites, which also integrates a number of protein signature databases as well as the appropriate search tools for them. The motifs found were then searched for in the InterPro database using InterProScan.

2D-structure determination

The distribution of the 2D structure in the motif region for each sequence was studied in order to propose a relation linking 2D structure to the function for each group of toxins. The 2D structure of all the sequences (training and test datasets) was determined based on the program PHD [36–38] using neural network approach, and the Softberry's software PSSfinder [39], which uses Markov chains probabilistic model.

Phylogenetic analysis

Phylogenetic analysis was carried out to:1) study the potential evolutionary relations between the new sequences and already known toxins, 2) infer potential common functions, and 3) classify the new ones (with unknown functions) with respect to the other toxins with known activities.

Since multiple-alignment is at the base of any phylogeny, alignments obtained from ClustalW results were used as entry for program PHYLIP [40] version 3.6 used to build phylogenetic trees. Distances were generated using Jones-Taylor-Thornton model. Phylogenetic trees were constructed using the neighbor-joining algorithm [41].

A construction of the phylogenetic trees for each of sequence group considered was also performed. Consensus trees were determined using the bootstrap method [42].

Statistical analysis

Correspondence analysis was performed, and biplot [43] was drawn, to illustrate the association between the motifs found and the activities of the toxins from the training dataset. To test the strength of this association, a χ² test of independence of factors was used where p-value was computed using Monte-Carlo simulation [44]. This procedure is implemented in the R software version 2.7 [45], and used in case where large sample theory is not valid (many cell having values less than 5). The same test χ² was used to measure association between the distribution of the 2D structure (determined by PHD and PSSfinder) of the motifs, and the toxins' activities. In order to arrange the toxins according to their activities (dependent variable) and the motifs determined and their 2D structures (structural variables), classification and regression trees (CART) [46] were built using the procedures implemented in the software Splus version 6.2 [47]. The built classification trees (for each 2D program results) were then used to predict the classes of the new sequences (test dataset), in a purely statistical way independent of the phylogeny results.

Motifs search results

The training dataset, in its different partitions, was used as input to GIBBS and MEME to conduct motif search. While GIBBS provides gapped or un-gapped motifs, MEME identifies only un-gapped motifs.

GIBBS results

MEME results

Using MEME and considering the relative diversity among K⁺ and Na⁺ channel effectors and the length range of these toxins, the maximum motif number was set to 6 (varied from 1 to 6) and maximum motif length was set to 20. Four motif runs were reasonable and provided the best and the most significant distribution of motifs among the groups of toxins. Therefore, we used this run as input to the MAST and Meta-MEME programs. These multilevel motifs are in PSSM format and the consensus (most probable) strings are shown in the Table 2 and Figure 1.

Motif Number	Group	Length	Consensus (Most Probable) Motif
1	Na+	20	GNACWCIELPDNVPIRIPGK
2	Na+	11	THLYEQAVVWP
3	birtoxin-like	20	NYPLDSSDDTYLCAPLGENP
4	K+	20	KDAGMRFGKCMNRKCHCTPK

Motif-1 was conserved in 32 sequences, all of them are α-type Na⁺ channel blockers. Motif-2 was conserved in 5 sequences, all of them are β-type Na⁺ channel blockers. Motif-3 was observed in 15 sequences among them 14 toxins belong to the birtoxin-like group, and only one is Na⁺ channel effector: CsvII. Motif-4 was observed in 16 sequences where 14 of them are active on the K⁺ channels and 2 are birtoxin-like peptides (birtoxin and ikitoxin). The remaining 14 sequences, (1 Na⁺ channel effector, and 13 K⁺ channel effectors) did not report statistically significant matches to any of the consensus motifs described in table 2. However, these results do not exclude that these sequences may contain other probabilistic variants of the motifs found. Indeed, 11 among the 13 K⁺ channel effectors were reported by GIBBS to have a conserved motif of length 20, that overlap with the one identified by MEME (motif4). The comparison of the results obtained by MEME with those obtained by GIBBS, showed that the motifs of length 20 determined by GIBBS were either included in the motifs determined by MEME or significantly overlapped with them. However, the results obtained by GIBBS are less significant in probabilistic terms and in the ability of differentiating among the groups of toxins. Therefore, only motifs obtained by MEME will be adopted for the rest of this work.

Identification of new K⁺ and Na⁺ channel effectors

Pattern and domain analysis

Pattern and domain analysis were further used to filter the results. All these 30 new K⁺ and Na⁺ channel effector sequences have at least one of the 4 conserved motifs previously found in this study. We aligned these 30 sequences and other known K⁺ and Na⁺ channel effectors. These motifs were overrepresented among all these sequences (newly identified and known toxins), thus, indicating that the identified motifs in the study are useful for searching potential new K⁺ and Na⁺ channel effectors from scorpions or plants.

InterPro results

Searching InterPro for the identified motifs showed that motif-1 corresponded to 2 hits: PD000908 characteristic of long chain scorpion toxins, and PF00537 corresponding to Toxin_3 domain characteristic of Na⁺ channel inhibitors from scorpion or plants. Motif-2 and motif-3 did not correspond to any match. Motif-4 corresponded to 2 hits: PD003586 characteristic of short chain scorpion toxins, and PF00451 corresponding to Toxin_2 domain. The fact that 2 motifs among the 4 identified overlap with other already reported in InterPro is a corroboration of our approach.

2D results

Using PHD and PSSfinder, all 2D structures of the sequences were determined. PHD results showed that the 2D structures displayed by motif-1, motif-2, motif-3, and motif-4 were mainly β-sheet, α-helix, β-sheet, and α-helix, respectively. Using PSSfinder, the 2D structure motif-1, motif-3, and motif-4 were mainly β-sheet, β-sheet, α-helix, or β-sheet, respectively. The 2D structure displayed by motif-2 was unidentified.

Phylogenetic analysis results

All 152 sequences (training and/or test datasets) including known and newly identified K⁺ and/or Na⁺ channel toxins, found in our analysis, were used to build phylogenetic consensus trees. Using the training dataset only (79 sequences), the consensus phylogenetic tree obtained revealed 3 major groups or clusters (figure 2a):

Group G1; contains all sequences of α-type Na⁺ channel toxins. Group G2 contains toxins having some activities on K⁺ channel. This group can be further subdivided into two subgroups: S1 containing toxins active on K⁺ channel only, and S2 involving birtoxin-like peptides. Group G3 contains β-type Na⁺ channel toxins. The group S2 of birtoxin-like peptides was situated between clusters S1 and G3, where S1 contains sequences active on K⁺ channel and G3 contains β-type Na⁺ channel toxins. Therefore, S2 can be considered as a transition group from toxins active on K⁺ channel to those β-type active on Na⁺ channel (figure 2b). This consensus tree presents the classical divergence between α and β Na⁺ channel toxins, and it is consistent with previous phylogenetic constructions made by [14, 49] even though, they did not use the same set of toxin sequences.

To infer the function of the test dataset toxins (72 sequences), all sequences from training and test datasets were used to build a consensus phylogenetic tree, see figure 2. We examined the whereabouts of the sequences from the test dataset with respect to the clusters (groups) defined by the training dataset. Thus, to perform a final check on the reliability of the results obtained, we examined the whereabouts of the 42 sequences of the test dataset with known function in the output of PHYLIP with respect to the clusters defined by the training dataset (79 sequences). The results show that these sequences (42) with known activities fall within compatible functionally defined clusters (clusters defined by toxins of the training dataset), thus, supporting our results.

The activities of the remaining 31 sequences were inferred in the same manner, see table 3 and figure 1. We can see that most newly identified peptides were grouped with sequences with known function toxins.

Statistical analysis results

The association between generated motifs and the toxins activities was very significant, simulated p-value < 0.0001. See biplot figure 3

Similarly, and while the majority of the motifs displayed a β-sheet structure, the association between motifs 2D structure distribution, and the protein classes was evaluated with the same χ² test, and simulated p-value < 0.001. Therefore, 2D structures of the motifs were incorporated as explanatory variables in the CART model.

Classification trees were adjusted using the training dataset of 79 sequences. Misclassification error rates for PHD and PSSfinder 2D based trees were less of 0.05, indicating good classifications. The same models were used to predict the activities of the 72 sequences identified by MAST. The 42 sequences with known activities were correctly classified, while the remaining 30 new sequences of the test dataset were assigned activities using CART model. The predicted activities were in perfect concordance with the results of the phylogeny. This fact confirms the classification tree models built using the training dataset as well as the results of the phylogeny.