Venom gland transcriptomics for identifying, cataloging, and characterizing venom proteins in snakes

Toxicon 93 (2015) 1e10

Contents lists available at ScienceDirect

Toxicon journal homepage: www.elsevier.com/locate/toxicon

Review

Venom gland transcriptomics for identifying, cataloging, and characterizing venom proteins in snakes Rajeev Kungur Brahma a, Ryan J.R. McCleary b, R. Manjunatha Kini b, c, d, Robin Doley a, * a

Department of Molecular Biology and Biotechnology, Tezpur University, Tezpur 784 028, Assam, India Department of Biological Sciences, Faculty of Science, National University of Singapore, Singapore, Singapore c Department of Biochemistry, Medical College of Virginia, Virginia Commonwealth University, Richmond, VA, USA d University of South Australia, School of Pharmacy and Medical Sciences, Adelaide, South Australia 5001, Australia b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 24 July 2014 Accepted 27 October 2014 Available online 29 October 2014

Snake venoms are cocktails of protein toxins that play important roles in capture and digestion of prey. Significant qualitative and quantitative variation in snake venom composition has been observed among and within species. Understanding these variations in protein components is instrumental in interpreting clinical symptoms during human envenomation and in searching for novel venom proteins with potential therapeutic applications. In the last decade, transcriptomic analyses of venom glands have helped in understanding the composition of various snake venoms in great detail. Here we review transcriptomic analysis as a powerful tool for understanding venom profile, variation and evolution. © 2014 Elsevier Ltd. All rights reserved.

Keywords: Snake venom Venom gland Transcriptome Venom complexity Venom protein cataloging Venom protein evolution Venomics

1. Introduction Venomous snakes are included in the Colubroidea (the so-called “advanced snakes”), and many possess well-developed venom delivery systems with tubular fangs and venom glands surrounded by muscle. Envenomation by these snakes can be a significant threat to human health (Chippaux and Theakston, 1987; Dolab et al., 2014; Gonzalez-Andrade and Chippaux, 2010; Gutierrez et al., 2006; Lalloo et al., 1996; Pierini et al., 1996; Ramos-Cerrillo et al., 2008; Vidal et al., 2007), especially those by front-fanged snakes such as members of the Viperidae and Elapidae, which are represented by over 600 species (Vonk et al., 2011). Some colubrid, lamprophiid, and natricid snakes also are venomous and, although many are not harmful to humans, some are potentially dangerous (Pyron et al., 2011, 2013). There also are potentially harmful species that are not known to cause human envenomation, either due to limits of morphology (small size, posteriorly-located fangs) or behavior (reluctance to bite). Traditionally, our knowledge of snake venoms was limited to species that cause human death and

* Corresponding author. E-mail address: [email protected] (R. Doley). http://dx.doi.org/10.1016/j.toxicon.2014.10.022 0041-0101/© 2014 Elsevier Ltd. All rights reserved.

debilitation, and even within those venoms, relatively few individual toxins have been fully examined. In recent years, significant progress has been made in understanding pathophysiological effects of envenomation and venom complexity by using cutting-edge proteomic and transcriptomic approaches. Such analyses have contributed in many ways including: 1) determining venom composition, 2) discerning variations in venom composition that inform our understanding of the evolution of snakes and their venoms, 3) identifying novel toxins that may help in the development of unique research tools and/or therapeutic agents, and 4) correlating specific pathological symptoms of envenomation with individual toxins. Snake venom is an evolutionary adaptation comprising many toxins that have evolved from cognate proteins (those which play a role in original, normal physiological function) through gene duplication followed by accelerated point mutation in protein coding regions (Nakashima et al., 1993; Ogawa et al., 1996; Ohno et al., 1998) and subsequent neofunctionalization (Brust et al., 2013; Fry, 2005; Fry and Wuster, 2004; Kini and Chan, 1999). The type of prey available is considered to be one of the main driving forces for evolution of snake venom composition (Daltry et al., 1996). Not only can there be variability among the venoms from snakes of separate genera, but also within individuals of the same

2

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

species (Chippaux et al., 1991). Factors indicated as influencing venom variation include age, sex (Menezes et al., 2006), locality (Durban et al., 2011; Jayanthi and Gowda, 1988; Minton and Weinstein, 1986), diet (Barlow et al., 2009; Creer et al., 2002; Daltry et al., 1996; Li et al., 2005b), season (Gubensek et al., 1974), frequency of extraction (Willemse et al., 1979) and even postextraction processing (Willemse and Hattingh, 1979). In fact, the venom from the same individual can be subject to ontogenetic change, seasonal variation, and differences in venom gland location (left vs. right) (Gubensek et al., 1974; Johnson et al., 1987). All these variations in venom compositions could be not only due to changes in the expression pattern of large numbers of isoforms within venom protein families (which may differ in their biological activities), but also to additional duplications, rapid ‘birth and death’ of toxin genes and accelerated evolution of exons (Fry et al., 2003b). Because of these factors, pathological symptoms in human victims may differ greatly even when envenomated by the same species. Emphasis traditionally has been placed on understanding either the most abundant or the most toxic proteins within snake venoms in order to assist in snake bite management, but the proteins found in smaller quantities were not considered due to difficulty in their purification and characterization. However, even such minor components may contribute to pathological symptoms either

directly or indirectly, through synergism. With recent advances in high-throughput techniques, the compositions of venoms can be deciphered with finer detail. Complementary DNA (cDNA) libraries of venom gland tissue and proteomic analyses of venoms are excellent tools for cataloging toxins including proteins of low abundance, thereby expanding our understanding of venom protein complexities (Dutertre et al., 2013; Francischetti et al., 2004; Jin et al., 2013; Lavergne et al., 2013; Siang et al., 2010). Subtle differences in venom protein structure and concomitant alterations in function (structureefunction relationships) need to be rigorously examined in the coming decades. With improvements in recombinant DNA technology, protein folding methods, high throughput screening, and structure determination, novel toxins can be characterized thoroughly and quickly for pharmacology and structureefunction relationships. These new capabilities have expanded our knowledge of the pharmacological profiles of venoms, helped produce better antivenoms through the field of antivenomics (Calvete et al., 2009; Calvete, 2010; Gutierrez et al., 2009; Harrison et al., 2007; Pla et al., 2012), aided in the identification of toxins that could be used as tools for understanding normal physiological processes, and developed as diagnostic tools, research means, or prototypes of therapeutic drugs. Further, such systems approaches help us decipher the evolutionary history, expression profiles and

Fig. 1. Strategies involved in transcriptomic analyses of snake venom glands. The red ‘X’ indicates cannot be used for cloning and expression. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

regulation of expression of venom toxin genes (Doley et al., 2009). The cDNAs of individual toxins or chimeras representing several toxins have also been used in raising antibodies that can neutralize their toxic effects (Azofeifa-Cordero et al., 2008). Thus, transcriptomics and proteomics have revolutionized the field of toxinology and opened several new avenues of research that were beyond our imagination only a decade ago. 2. Venom gland transcriptome generation 2.1. Based on cloning technology For successful construction of a venom gland cDNA library, it is imperative to have high levels of mRNAs encoding toxin genes. It has been well established that the glandular epithelial cells become shorter and quiescent when the glands are full, but become enlarged and highly active in the production of toxins when the venom glands are empty (De Lucca et al., 1974). After venom extraction, an empty venom gland requires 40e50 days to completely refill, but maximal mRNA concentrations have been observed at 4e8 days post extraction (Kerchove et al., 2004). Therefore, venom glands should be harvested 4e5 days after venom extraction, the time at which toxin mRNAs are being transcribed at the highest rate (Currier et al., 2012). Immediately after dissection of the venom glands, they are stored in an RNA stabilization reagent (such as RNAlater, Invitrogen) to avoid RNA (specifically mRNA) degradation. Venom gland tissue is homogenized, and the total RNA is isolated using one of many different methods. The integrity of the total RNA isolated can be determined by the presence of bright bands of ribosomal 28S and 18S RNA at approximately 4.5 and 1.9 kb, respectively, on a 1% agarose gel (Fig. 1). A sample can also be run through a Bioanalyzer (Agilent Technologies) to give quantitative and qualitative (RINdRNA integrity number) assessments of the sample (Fig. 1). mRNAs are isolated from total RNA utilizing a primer containing a polyT (15e20 bp) section with a 30 VN (V being any base other than T and N being any base) and a 50 known adaptor sequence. After this primer anneals, a reverse transcriptase (RT) is utilized to create a DNA strand complementary to the RNA and to add an overhang of a few additional bases. A second adaptor sequence will anneal to this overhang, and the RT will continue creating the complementary strand. The result is a single stranded cDNA with adaptor sequences on both ends that can be used for creation of

3

double stranded cDNA and later PCR amplification. Most of the commercially available cDNA construction kits provide vectors that contain restriction sites for ease of ligation of the gene of interest. Recombinant vectors containing the cDNA inserts are transformed into competent Escherichia coli cells and plated on LB agar plates containing ampicillin, IPTG, and X-gal for blue/white screening. Individual colonies represent one cDNA sequence, and the presence of an insert is confirmed by restriction digestion. Plasmids containing the cDNAs are isolated and subjected to sequencing (such as by the BigDye terminator reaction, Life Technologies) using universal primers. Low quality sequences and short expressed sequence tags (ESTs < 150 bp) can be removed using the Phred program (www. phrap.com) and adapter and vector sequences using the NCBI VecScreen (http://www.ncbi.nlm.nih.gov/) or CrossMatch programs. The ESTs obtained after removing the vector sequences can be assembled into contiguous clusters (contigs) using the CAP3 program (Huang and Madan, 1999). All contigs and singletons are subjected to BLAST searches against the non-redundant database of NCBI (http://www.ncbi.nlm.nih.gov/) using BLASTX and BLASTN algorithms (e-value cutoff < 105 with a minimum coverage of 100 bp and >90% identity) for identification of putative genes (Pahari et al., 2007b). The sequences that are not identified by BLAST or are of unknown function can be analyzed for presence of signal peptides using SignalP (http://www.cbs.dtu.dk/). This allows identification of the putative mature protein for further analysis and indicates whether it is a secretory protein. An annotation table containing all the relevant information about contigs and singlets can be generated followed by submission to GenBank dbEST or any other database. However, the criteria for manual analysis of the sequence is based on 1) similarity of transcripts to known toxin sequences available in the database, 2) dissimilar open reading frames to those found in databases, 3) presence of a signal peptide (as venom proteins are secretory proteins) and cysteine residues (most venom toxins are cysteine rich). Moreover, the presence of signature motifs unique to certain snake venom protein families can be useful for manual cataloging of the transcripts (Table 1). Further sequence analysis needs to be completed to check for nucleotide substitutions, deletions, and/or additions in coding regions, as many toxin families are rich in isoenzymes with differences in amino acid sequence.

Table 1 Some snake venom protein families and their signature motifs. Protein family abbreviations are as in the text. An “X” in the signature motif indicates any amino acid residue. Protein family

Approximate molecular weight range (kDa)

Signature motif

References

3FTx (long chain) 3FTx (short chain) 3FTx (non-conventional) PLA2 (group I)

6e9 6e9 6e9 13e15

(Menez, 1998; Tsetlin, 1999) (Menez, 1998; Tsetlin, 1999) (Menez, 1998; Tsetlin, 1999) (Scott and Sigler, 1994)

PLA2 (group II)

13e15

Kunitz-type serine protease inhibitor Snaclec

~7

Disulphide bonds: C1eC3, C2eC6, C4eC5, C7eC8 and C9eC10 Disulphide bonds: C1eC3, C2eC4, C5eC6, and C7eC8 Disulphide bonds: C1eC5, C2eC3, C4eC6, C7eC8 and C9eC10 Conserved cysteines and bonds: C1eC8, C2eC14, C3eC5, C4eC13, C6eC12, C7eC10, and C9eC11 Disulphide bonds: C1eC13, C2eC4, C3eC12, C5eC14, C6eC11, C7eC9, and C8eC10 Conserved cysteines: C[8X]C[15X]C[4X]YGGC[12X]C[3X]C

SV-LAAO SVMP SVSP CRISP Hyaluronidase

26e28 57e68 (monomer) 110e125 (dimer) 36e60 36 21e29 33e110

(Scott and Sigler, 1994) (Gojobori and Ikeo, 1994)

Carbohydrate recognition domain (CRD): C31, G69, T92, P97, C106, A120, C123 and C131 FAD binding domain, substrate binding domain and helical domain

(Hirabayashi et al., 1991) (Pawelek et al., 2000)

Zinc binding domain Catalytic triad: S195, H57 and D102 Cysteine rich C-terminal domain Conserved cysteines: C340, C211, C227, C365, C370, C376 and C429

(Fox and Serrano, 2005) (Serrano and Maroun, 2005) (Guo et al., 2005) (Harrison et al., 2007)

4

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

Table 2 List of snake species (and corresponding family) for which cDNA libraries have been compiled. This list does not include studies that utilized transcriptomic techniques to analyze only specific toxins or families of toxins. Platform refers to the method by which the cDNA library was compiled, whether through clonal (cDNA library) or through next-generation sequencing (Illumina and 454) techniques. Taxonomy is per Pyron et al. (2011) and reptile-database.org. Family and Species Colubridae Dispholidus typus Telescopus dhara Thrasops jacksonii Trimorphodon biscutatus Dipsadidae Liophis (¼Erythrolamprus) miliaris Liophis (¼Erythrolamprus) poecilogyrus Philodryas olfersii

Thamnodynastes strigatus Elapidae Austrelaps labialis Bungarus flaviceps Bungarus multicinctus Demansia vestigiata Drysdalia coronoides Micrurus altirostris Micrurus corallinus Micrurus fulvius Naja atra Oxyuranus microlepidotus Homalopsidae Cerberus rynchops Enhydris polylepis Lamprophiidae Leioheterodon madagascariensis Psammophis mossambicus Natricidae Rhabdophis tigrinus Viperidae Agkistrodon piscivorus leucostoma Atropoides mexicanus Atropoides picadoi Azemiops feae Bitis gabonica Bothriechis lateralis Bothriechis schlegelii Bothrops alternatus Bothrops asper Bothrops atrox Bothrops insularis

Bothrops jararaca

Platform

Reference

cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Fry et al., 2008) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Fry et al., 2012) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Ching et al., 2006; Fry et al., 2008; Fry et al., 2012) cDNA library (Ching et al., 2012) cDNA cDNA cDNA cDNA cDNA cDNA cDNA

library library library library library library library

(Doley et al., 2008b) (Siang et al., 2010) (Jiang et al., 2011) (St Pierre et al., 2007) (Chatrath et al., 2011) (Correa-Netto et al., 2011) (Correa-Netto et al., 2011; Leao et al., 2009) Illumina (Margres et al., 2013) cDNA library (Jiang et al., 2011) cDNA library (Fry et al., 2008) cDNA library (OmPraba et al., 2010) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Jia et al., 2008) 454 (Durban et al., 2011) 454 (Durban et al., 2011) cDNA library (Fry et al., 2008; Fry et al., 2012) cDNA library (Francischetti et al., 2004) 454 (Durban et al., 2011) 454 (Durban et al., 2011) cDNA library (Cardoso et al., 2010) 454 (Durban et al., 2011) cDNA library (Neiva et al., 2009) cDNA library (Junqueira-de-Azevedo and Ho, 2002; Valente et al., 2009) cDNA library (Cidade et al., 2006; Zelanis et al., 2012) cDNA library (Rodrigues et al., 2012)

Bothropoides (¼Bothrops) pauloensis Bothrops jararacussu Causus rhombeatus Cerrophidion godmani Crotalus adamanteus

cDNA library cDNA library 454 454, Illumina

Crotalus durissus collilineatus (¼terrificus) Crotalus horridus

Illumina

(Kashima et al., 2004) (Fry et al., 2008) (Durban et al., 2011) (Rokyta et al., 2011; Rokyta et al., 2012) cDNA library (Boldrini-Franca et al., 2009) (Rokyta et al., 2013)

Table 2 (continued ) Family and Species

Platform

Reference

Crotalus simus

454

Deinagkistrodon acutus Echis carinatus shochureki Echis coloratus Echis ocellatus

cDNA cDNA cDNA cDNA

(Durban et al., 2011; Durban et al., 2013) (Zhang et al., 2006) (Casewell et al., 2009) (Casewell et al., 2009) (Wagstaff et al., 2009; Wagstaff and Harrison, 2006) (Casewell et al., 2009) (Junqueira-de-Azevedo et al., 2006) (Pahari et al., 2007b)

library library library library

Echis pyramidum leakeyi cDNA library Lachesis muta cDNA library Sistrurus catenatus edwardsii cDNA library

2.2. Based on next generation sequencing (NGS) technologies Early transcriptomic methods, such as candidate gene-based Northern blot analysis (Alwine et al., 1977) and Sanger sequencing (Sanger et al., 1977) were slow and expensive, while new technologies have increased throughput and reduced the cost. Each of these powerful technologies has benefits and limitations that make it important to know the desired data outcome before using them (Table 2). For high-throughput sequencing of transcripts using NGS technologies, the sample must be prepared with utmost care, since cross-contamination is difficult to detect. Tissue collection and RNA extraction are basically as described (Section 2.1). For sequencing, the initial RNA or resulting dsDNA must be fragmented, which is generally achieved by utilizing random primers (for RNA prior to cDNA creation), nebulization (ds cDNA) or enzymatic digestion with highly diluted DNaseI in the presence of Mg2þ (ds cDNA). The DNA fragment ends are polished/repaired (i.e., protruding 30 and 50 ends are removed or filled in), and platform-specific upstream and downstream adapters are ligated to both ends. The adapter ligated DNA fragments are size selected (500e800 bp) (Rokyta et al., 2011; Schwartz et al., 2010) on agarose gels. The library fragments are annealed or impregnated on capture beads (FLX and GA systems) or attached to oligo-derivatized surfaces of flow cells (SOLiD system) and prepared for amplification. For high-throughput sequencing, the 454 system uses pyrosequencing (chemiluminescence-based detection of pyrophosphates released during incorporation of nucleotides), the Ion Torrent uses semi-conductor sequencing (detection of hydrogen ions released during sequencing), and the SOLiD and 454 systems utilizes emulsion PCR ligation methods (Droege and Hill, 2008; Mardis, 2008). The latter two use emulsion droplets to isolate single DNA templates in separate micro reactors where amplification is accomplished. The Illumina platform uses bridge amplification, in which DNA molecules attached to a solid surface are amplified in situ, generating clusters of identical DNA copies (Bentley, 2006; Mardis, 2008). Two other current platforms (Pacific Biosystems and Helicos Biosystems) do not utilize amplification of template DNA for sequencing, but rather use single molecules of DNA (Ozsolak et al., 2009). In the case of Pacific Biosystems, this yields potential read lengths of over 10,000 bp (Eid et al., 2009). NGS technologies have aided in the construction of snake venom transcriptomes with greater coverage in shorter periods of time (Wang et al., 2009). The largest venom gland transcriptome assembled by NGS thus far is for Deinagkistrodon acutus (Chinese moccasin) with 8696 ESTs (Zhang et al., 2006). Rokyta and colleagues (Rokyta et al., 2011) generated a venom gland transcriptome of the Eastern diamondback rattlesnake (Crotalus adamanteus) using Roche 454 sequencing technology. Durban and colleagues (Durban et al., 2011) investigated the venom gland transcriptomes of eight Costa Rican snakes e Bothrops asper (terciopelo; from Caribbean and Pacific populations), Bothriechis

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

lateralis (side-striped palm viper), Bothriechis schlegelii (eyelash viper), Atropoides picadoi (Picado's pit viper), Atropoides mexicanus (jumping viper), Crotalus simus (Central American rattlesnake), and Cerrophidion godmani (Godman's montane pit viper) e using the 454 Genome Sequencer FLX system. They simultaneously sequenced eight snake venom gland transcriptomes in two runs, demonstrating large scale transcriptome coverage within a short time using NGS technologies (Durban et al., 2011). Table 2 is a list of current snake venom gland transcriptomes created using both general techniques.

3. A brief history of venom gland transcriptomics Early on in attempts to categorize toxin genes, it was noted that polyadenylated mRNAs could be isolated from venom gland tissue (Vandenplas et al., 1985). This was an important step because it indicated the potential for determining toxin gene sequences. Although this study yielded only a few mRNA species, the authors indicated that their techniques may not be sensitive enough to detect all mRNAs present. In the same year, the cDNA sequence of erabutoxin a (a three finger toxin e 3FTx) from Laticauda semifasciata was publisheddthe first such sequence from a snake venom gland mRNA (Tamiya et al., 1985). Subsequently, a number of individual toxin gene sequences were reported including a phospholipase A2 (PLA2) from Aipysurus laevis (Ducancel et al., 1988); a PLA2 from Notechis scutatus scutatus (Ducancel et al., 1988); crotoxin, a heterodimeric PLA2 complex from Crotalus durissus terrificus (Bouchier et al., 1988); and ammodytoxin, a PLA2 from Vipera ammodytes (Pungercar et al., 1991); followed by the first attempts at global cataloging of venom gland transcripts. The first venom gland cDNA library was compiled using cloning technology and was for the viperid Bothrops insularis (Junqueirade-Azevedo and Ho, 2002). The resulting cDNA library yielded 345 clones with toxin related sequences that represented 85 different clusters (unique sequences). These included PLA2s, snake venom C-type lectin-like proteins (Snaclecs), snake venom L-amino acid oxidases (SV-LAAOs), snake venom metalloproteases (SVMPs), snake venom serine proteases (SVSPs), a cysteine-rich secretory protein (CRISPs, formerly known as helveprins), bradykinin potentiating peptide (BPP) precursors, a vascular endothelial growth factor, and a nerve growth factor. The 85 toxin-related clusters represented 28.6% of the unique sequences, and the toxin genes in total represented 56.6% of all resulting sequences. This showed that (a) toxin genes were a large percentage of expressed genes in the venom gland, (b) the potential diversity of toxin genes was higher than previously known, and (c) such techniques can lead to the discovery of novel toxin genes and families. Since the first clonally-produced cDNA library, many others have been produced using both traditional cloning techniques and newer powerful NGS platforms (see Table 2). Cloning-based transcriptomes have been produced even recently (Rodrigues et al., 2012) for Bothropoides pauloensis, and this technique is still a viable alternative to NGS, because it is relatively inexpensive to produce 1000 to 2000 sequences and only requires equipment and reagents readily available in molecular biology laboratories. The first use of NGS for cataloging toxin genes from venom glands was done fairly recently using 454 technology to catalog toxins from the eastern diamondback rattlesnake (C. adamanteus) (Rokyta et al., 2011). They found 40 toxin sequences from 13 different toxin families, including two that had not previously been seen in crotalid snakeda putative phospholipase B and a putative epidermal growth factor. This study further exhibited the ability to detect toxins in very low abundance in the mRNA pool, something that is more difficult to do using cloning techniques.

5

4. Transcriptomics for cataloging toxin genes and potential venom proteins Snake venoms can contain toxins belonging to many different protein families, some examples of which include the PLA2s, 3FTxs, Kunitz-type serine protease inhibitors (Kunitz), snaclecs, SV-LAAOs, SVMPs, SVSPs, and CRISPs. Although all venomous snakes belong to the superfamily Colubroidea, some toxin families are more prevalent in certain taxa than in others. Cataloging of toxin genes via venom gland transcriptomics (as discussed earlier) helps in predicting the presence of toxins from various protein families in the venom. Such species-specific toxin profiles help in understanding pharmacological properties and pathological symptoms. They also allow for comparisons ranging from the intra-individual to interfamily level and give a greater understanding of the diversity of toxins in nature. With more venom profile data available, it is easier to predict the presence of related toxins using phylogeny-based approaches. Toxin cataloging can be useful in either targeted searches for individual toxins or isoforms within one toxin family or in attempts to discover all toxins within a species (global cataloging). For the former, specific primer sequences designed based on previously published sequences (e.g., http://www.ncbi.nlm.nih.gov/nuccore) can be utilized to amplify genomic DNA (gDNA) or complementary DNA (cDNA) constructs. In general, the signal peptide and 30 untranslated sequences of toxin genes are highly conserved across snake taxa. Simple PCR followed by agarose gel electrophoresis indicates the presence of a type of toxin, although it may not indicate the diversity of such toxins actually present. These genes can subsequently be cloned and sequenced for detailed comparative studies. Attempts at discovering the entire inventory of toxin genes within a venom gland also have traditionally utilized cloning techniques. However, the advent of NGS technologies has dramatically increased the efficiency of sequencing. While a cDNA library based on cloning generally includes up to a few thousand sequences, NGS can yield hundreds of thousands to millions of sequences. Although NGS drastically reduces the time and effort needed to generate data from tissue, it creates the need for competency in bioinformatics in order to logically and logistically analyze such data to achieve meaningful results. 5. Transcriptomics for studying toxinology 5.1. Documenting toxin profiles in venoms Probably the most powerful use of transcriptomics related to venomous animals is in their use to examine the relative types and amounts of toxin genes present in any given venom gland. This profiling allows for the discovery of novel toxins and toxin families as well as for comparisons with other species. Whether done via traditional cloning or more recent NGS techniques, the process yields useful information for further characterization of interesting new toxin genes or for understanding the relationships among species in a venom-oriented context. As an early yet fairly recent example of profiling, the transcriptome of the South American Bothrops jararacussu was created utilizing cloning techniques and yielded 549 ESTs, of which 340 (62%) were toxin sequences (Kashima et al., 2004). Overall, nine different toxin families were represented, with PLA2s being 197 of the toxin ESTs. Interestingly, this study also examined the relative presence of three different PLA2 sub-types (basic Lys49e83.2%; basic Asp49e8.1%; and acidic Asp49e8.6%). A second example published at nearly the same time was for the African Bitis gabonica, another viperid snake (Francischetti et al., 2004). In this cDNA library, 12 toxin families were represented utilizing sequences

6

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

combined into 300 unique contigs, of which the major toxin family was the SVMPs (30%). A simple comparison of the overall percentages of these two viperid species shows major differences in overall toxin composition, which is extremely useful for understanding variation among and within different taxonomic levels. More recently, NGS technologies have been employed for toxin profiling. Using the 454 platform for snake venom gland transcriptomics, a comparative study of seven species of viperid snakes from five different genera were characterized (Durban et al., 2011). The results indicate some major taxonomic differences within a family that is considered strongly monophyletic. Differences ranged from minor to major alterations in relative percentages of toxin families represented. Further, because NGS techniques allow for discovery of rare transcripts, there were many cases of presence/ absence among the taxa represented. A closer examination of two geographically-distinct B. asper individuals showed relatively similar toxin family patterns with the major exception that the ohanins were present as 2.8% of the protein family hits for the individual from the Caribbean region, but 11.9% of the hits for the individual from the Pacific region. Such studies show that transcriptomics can give a good starting point for characterization, but also have much potential to uncover the natural variation of toxins from snake venom.

treatment. It has been known for some time that ontogenetic shifts in venom composition occur in snakes (Gutierrez et al., 1990). However, only recently this phenomenon has been examined at the molecular level. A cDNA library comparison between neonate and adult Bothrops jararaca using cloning techniques showed a reduction in the relative percentage of SVMP transcripts from neonate (53.2%) to adult (29.9%). This was accompanied by increases in BPP precursors (9.6e23.2%), snaclecs (10.8e22.3%), and SVSPs (3.3e8.1%) (Zelanis et al., 2012). These temporal expression patterns reflect venom transcriptome complexity and potential shifts in the pharmacological activities of the resulting secreted venom. It is known that the complexity of snake venom within a single species can be extremely high. However, the real extent of this is not known for most species of snakes. As an example, analysis of the venom gland transcriptome of E. ocellatus (African saw-scaled viper), a medically important snake in West Africa, revealed a diverse expression of SVMPs and also several distinct isoforms of BPPs, PLA2s, Snaclecs, SVSPs and SV-LAAOs (Wagstaff and Harrison, 2006). This study uncovered a venom more complex than previously thought for this species, which is important because the presence of a diverse group of SVMPs has been indicated as causing the severe hemorrhagic conditions observed in humans envenomated by E. ocellatus, and coagulopathy might be due to the SVSPs and Snaclecs found in the venom (Wagstaff and Harrison, 2006).

5.2. Identifying novel toxins 5.4. Developing suitable antivenoms Transcriptome analyses are useful in identifying novel toxins in existing families as well as new families of toxins. 3FTxs were previously known to be present only in the venoms of elapids (elapinae-cobras, kraits, and mambas; hydrophiinae-sea snakes), and colubrids, however their presence in the viperid venom was first reported through trancriptomic analysis of Lachesis muta (Junqueira-de-Azevedo et al., 2006). Similarly transcripts belonging to non-conventional 3FTxs were also identified in the venom gland transcriptome of Sistrurus catenatus edwardsii (Pahari et al., 2007b). A new class of 3FTx was found through a partial cDNA library of Ophiophagus hannah (king cobra) venom gland tissue. This protein has been reported to act as a natural exogenous beta-blocker (Rajagopalan et al., 2007). Phospholipase B and renin-like aspartic proteases, both new venom toxin families, were reported for the first time in the venom glands of Drysdalia coronoides (Chatrath et al., 2011) and Echis ocellatus (Wagstaff and Harrison, 2006), respectively. Kazal-type inhibitor-like proteins were reported in the transcriptome of B. schlegelii and B. lateralis (Durban et al., 2011). Further transcriptomics can be used in studying venom profiles of snakes that produce small quantities of low-toxicity venom as well as in understanding the diversity of genes within a snake family (Casewell et al., 2009). For example, we examined the transcriptome of the venom gland of Cerberus rynchops, a homalopsid snake, using clonal techniques (OmPraba et al., 2010). This study found only transcripts for members of four venom protein families: SVMPs (30%), snaclecs (22%), CRISPs (22%), and the venom ficolinlike proteins (veficolins, 26%). Further comparison with proteomic data (via tandem mass spectrometry after tryptic digest of crude venom) showed a similar pattern. This was the first study to identify the veficolins, indicating that even species with low transcriptomic complexity may produce novel toxins worthy of further study. 5.3. Cataloging toxin variations Transcriptomic cataloging of toxin genes expressed in the venom gland of a given snake elucidates natural variation in venom components and potential pharmacological effects of snake bite, and thus assists in developing suitable strategies for snake bite

An important practical use of venom gland transcriptomics is in development of antivenom for snakes of medical importance. The South American painted coral snake (Micrurus corallinus) has venom that is highly neurotoxic and known to cause deaths in humans due to diaphragm paralysis. Although bites by this species are fairly rare, they are potentially severe, and antivenom is difficult to obtain due to low venom yields and difficulty in maintaining this animal in captivity. In order to understand venom complexity and to identify potential antigens for antivenom development, a venom gland cDNA library was constructed from 1438 ESTs (Leao et al., 2009). This library included a low diversity of pre-synaptic PLA2s and an extremely high diversity of post-synaptic 3FTxs, both at seemingly high expression levels. Further, the transcriptome identified eight other toxin families, indicating a diverse venom profile that can potentially be used for developing antivenom. 5.5. Studying venom evolution cDNA sequence analysis can give great insight into not only what toxin genes are present (and, therefore, what proteins may be present in the venom), but also their evolutionary history. Such analysis can help to elucidate certain relationships among different species and can also help in understanding how certain toxins came to be present or absent in the venom. The following transcriptomic examples have given insight into some of the evolutionary processes in action. Transcriptomes can give insight into the loss of venom toxicity, a phenomenon likely to have happened many times in colubroid snakes. Sea snakes (Hydrophiinae) generally are known to have extremely toxic venom, and the thought is that such strong venom is needed to rapidly subdue fish prey so that they cannot swim away from the snake before dying (Li et al., 2005a, 2005b). In a marine environment, a slowly dying fish could swim to an inaccessible location or be consumed by another nearby predator, thereby wasting the snake's venom. However, at least one species, Aipysurus eydouxii, is known to feed exclusively on fish eggs (Voris and Voris, 1983), basically negating the need for venom. As venom is thought to be energetically costly to produce (McCue, 2006), its

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

presence in A. eydouxii is counterintuitive. The cDNA library of A. eydouxii includes 31% housekeeping genes, 42% unknowns, and 27% toxin genes (21% PLA2s, 5% 3FTxs and 1% toxin-like proteins) (Li et al., 2005a). It was observed that the sequence of the only 3FTx has a dinucleotide deletion that leads to premature truncation of the toxin and a loss of functional neurotoxicity. Because the dominant snake venom protein families in the sea snakes are 3FTxs and PLA2s (Fry et al., 2003a), the functional loss of this 3FTx in A. eydouxii venom leads to a concomitant 40e100 fold decrease in venom toxicity (Tu, 1973). Thus transcriptomics can help detect potential mechanisms by which venom toxicity is lost. Pahari et al. (Pahari et al., 2007a) showed that sea snake venoms have the lowest, sea kraits (Elapinae, formerly Laticaudinae) have greater, and terrestrial elapids (Elapinae) have the highest relative venom complexity, which could be related to their respective prey diversity. Both 3FTx and PLA2 in the two sea snakes, Lapemis curtus (generalist feeder) and Acalyptophis peronii (restricted feeder), are simple and do not possess notable diversity, suggesting that the venoms genes are well conserved despite the large adaptive radiation of the two snakes. Ecological variations might have played an important role in determining the evolutionary trends of the two snakes along ecological niches (specialist and generalist) and not completely because of a distant phyletic relationship between them. Expression profile of these snakes suggests that the evolutionary trends of toxins may be linked to their diet (Pahari et al., 2007a). In a similar study of the Australian copperheads, represented by three species e Austrelaps superbus (lowland copperhead), Austrelaps ramsayi (highland copperhead), and Austrelaps labialis (pygmy copperhead), transcriptomic analysis indicated evolutionarily divergence in venom composition. A. labialis is smaller in size than the other species and mainly feeds on small animals (Shine, 1987). The cDNA library of the A. labialis venom gland includes genes for six venom protein families: 3FTx (45%), PLA2 (33%), Kunitz (9%), CRISP (8%), SVMP (3%) and snaclec (2%) (Doley et al., 2008b). Interestingly, 40% of the toxin gene transcripts are truncated due to addition or deletion of nucleotides. This likely leads to lower toxicity of the venom, and the combination of difference in toxicity and smaller size may have resulted from the partitioning or specialization of prey types in the ancestral form of the snake. It also has been reported that, out of the 19 neurotoxin isoforms observed in the venom gland cDNA library, six have prematurely truncated transcripts (Doley et al., 2008b). Toxin genes in the venom gland are thought to originate initially through gene duplication that allows for reduced selective pressure on one isoform of the gene that can then mutate, yielding a protein with diverse structure and function (Chang et al., 1999; Kordis and Gubensek, 2000; Moura-da-Silva et al., 1995). This novel gene can then be upregulated in the venom gland and secreted as a potential toxin. Transcriptomics methods have been utilized to examine the first step of this evolution. Based on sequence homology, group C prothrombin activators from elapid venoms are similar to factors Xa and Va of the blood coagulation cascade and are likely to have originated through gene duplication events (Rao et al., 2004). A comparative study of prothrombin activator genes from the Pseudonaja textilis venom gland and normally-functioning factor X from its liver yielded an apparent intermediary gene (Reza et al., 2006). This gene, whose primary sequence is similar to pseutarin C, is not expressed in the venom (indeed, it is only expressed in low amounts in the liver). Recently, based on the cDNA library of S. c. edwardsii, Doley et al. deduced the gene structure of a novel protein called ku-wap-fusin (Kunitz waprin fusion protein) (Doley et al., 2010). After analyzing the cDNA sequence and gene structure, they hypothesized that this fused protein might be an evolutionary intermediate of the Kunitz

7

toxins and the waprins. These two previous examples indicate intermediate evolutionary steps between gene duplication and eventual recruitment as a secreted protein. Similar cDNA sequences were also found in Sistrurus catenatus tergeminus venom gland (Unpublished observations). Such evolutionary mechanisms might be contributing to venom complexity and, therefore, toxicity. b-bungarotoxin is the principle toxic component of the venom of many Bungarus species, and it is a pre-synaptic toxin. It is a covalent heterodimeric complex between a PLA2 and a Kunitz (MacDermot et al., 1978). Transcriptomic analysis of the venom gland of Bungarus flaviceps showed that 2.87% of the transcripts were of the A chain and 32.01% were of the B chain of b-bungarotoxin. This suggests that b-bungarotoxin is likely to be a major protein in the venom, and it is known to be the main lethal component (Siang et al., 2010). In addition, a cluster of Kunitz toxins similar to the B chain of b-bungarotoxin, but without the extra cysteine required for covalent interaction with the A chain, was also observed (Siang et al., 2010). As the cysteine residue (when present) is located on the C-terminal end of the mature protein, the addition or loss of such a residue can help lead to or negate (respectively) the formation of a dimer. Cysteine residues in toxins cause internal disulfide bonds and are thought to lead to their stability. Toxins with odd numbers of cysteine residues, therefore, have an unpaired thiol group that can potentially bond with another protein, thereby creating a dimer. Irditoxin, reported from Boiga irregularis (brown tree snake), is a novel covalently-linked heterodimeric 3FTx (Pawlak et al., 2009). Analysis of the cDNA sequence of both its subunits reveals that the extra cysteine in each subunit is due to substitution of nucleotides. These extra cysteine residues allow for the formation of a heterodimeric complex that specifically targets the avian neuromuscular junction. A further potential mechanism for increasing complexity in venoms has also been examined with transcriptomics. Using the cDNA of venom gland tissue from S. c. edwardsii (desert massasauga), gene structures of viperid 3FTxs were reported for the first time (Doley et al., 2008a). Based on gene sequences, viperid 3FTxs appear to have evolved through a phenomenon called accelerated segment switch in exons to alter targeting (ASSET). It has been hypothesized that some multi-nucleotide sequences (segments) undergo switching, leading to changes in the amino acid sequence that alters the targeting of various ion-channels/receptors. This phenomenon, initially discovered through transcriptomic analysis, has also been seen in other taxa (Doley et al., 2009). Recently Sunagar et al., have proposed a new hypothesis to explain the molecular evolution of the 3FTx gene called “Rapid Accumulation of Variations in Exposed Residues” (RAVER) by analyzing the nucleotide sequences (Sunagar et al., 2013). Accordingly, hypermutable sites occur on surfaces of 3FTxs, which undergo changes without disrupting the structural integrity. Amino acids of this region undergo rapid mutation and are selected through positive selection if functionally relevant or else get eliminated through negative selection (Sunagar et al., 2013). Though many theories have been put forward to understand the molecular mechanism of snake venom proteins still the puzzle remains to be solved. With more transcriptomic and genomic understanding of the venom proteins, the molecular mechanism of venom protein evolution will be uncovered. 6. Conclusion and future directions Transcriptomic profiling is an excellent tool for understanding the composition and complexity of snake venoms. With recent advances in high throughput technologies, these complexities have been revealed in greater detail. Such studies have contributed to

8

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

predicting the pharmacological profiles of snake venoms that otherwise require proteomic approaches and animal studies. Variations in venom composition and relative expression of venom proteins can be elucidated using this technology. Moreover, mechanisms of venom protein evolution and neofunctionalization, which further contributes to venom complexity, have been deciphered. Identification of previously ignored proteins of low abundance has increased rapidly and could contribute to understanding venom toxicity and identifying potential therapeutic leads in the near future. Identification of toxin gene structures has helped in understanding potential function and in dissecting structureefunction relationships. However, such studies do not provide information on posttranslational modifications such as glycosylation and proteolytic cleavage. They also fail to provide information on dimerization or indicate exactly which mRNA transcripts are likely to be translated into protein. Moreover, currently animals need to be sacrificed for extraction of the venom gland tissue (although recent studies have shown that lyophilized venom contains at least some useful RNA) (Currier et al., 2012). Hence transcriptomic analysis needs to be complimented with proteomic analysis for full pharmacological profiling of snake venoms (OmPraba et al., 2010). However, by using transcriptomic analysis the limitations of antivenom therapy, which is currently the only available effective treatment for envenomation, can be overcome. Utilizing high-throughput transcriptomic analyses, regional variations in venom composition can be identified, which will lead to design and improvement of specific and efficacious anti-venoms. Further, unique toxins in the venoms of snakes from particular geographic locations can be identified and exploited for preparation of immunological detection kits. This is extremely important, as in most instances the snake species responsible for any given bite remains unidentified. Although some of the typical processes involved in creating cDNA libraries and transcriptomes have been included here, the technology has evolved rapidly over the past decade. This has been propelled by the need for high quality genetic data, and has in turn led to innovations in bioinformatics analysis of incredibly large datasets. As such large datasets become the standard for transcriptomic studies, the technology will likely evolve further, decreasing the time and cost involved in producing and analyzing such data. There have already been a few iterations of this cycle, and it is expected to continue further, which will lead to asking more complex questions in the future. Ethical statement This review does not involve any animal experiments. Conflict of interest We declare that there is no conflict of interest. Acknowledgments RKB acknowledges Department of Biotechnology, Govt. of India for the Junior Research Fellowship. The authors acknowledge Dr. Rupak Mukhopadhyay for his critical comments and suggestions. References Alwine, J.C., Kemp, D.J., et al., 1977. Method for detection of specific RNAs in agarose gels by transfer to diazobenzyloxymethyl-paper and hybridization with DNA probes. Proc. Natl. Acad. Sci. U. S. A. 74, 5350e5354. Azofeifa-Cordero, G., Arce-Estrada, V., et al., 2008. Immunization with cDNA of a novel P-III type metalloproteinase from the rattlesnake Crotalus durissus

durissus elicits antibodies which neutralize 69% of the hemorrhage induced by the whole venom. Toxicon 52, 302e308. Barlow, A., Pook, C.E., et al., 2009. Coevolution of diet and prey-specific venom activity supports the role of selection in snake venom evolution. Proc. Biol. Sci. 276, 2443e2449. Bentley, D.R., 2006. Whole-genome re-sequencing. Curr. Opin. Genet. Dev. 16, 545e552. Boldrini-Franca, J., Rodrigues, R.S., et al., 2009. Crotalus durissus collilineatus venom gland transcriptome: analysis of gene expression profile. Biochimie 91, 586e595. Bouchier, C., Ducancel, F., et al., 1988. Cloning and sequencing of cDNAs encoding the two subunits of Crotoxin. Nucleic Acids Res. 16, 9050. Brust, A., Sunagar, K., et al., 2013. Differential evolution and neofunctionalization of snake venom metalloprotease domains. Mol. Cell. Proteomics 12, 651e663. Calvete, J.J., 2010. Antivenomics and venom phenotyping: a marriage of convenience to address the performance and range of clinical use of antivenoms. Toxicon 56, 1284e1291. Calvete, J.J., Sanz, L., et al., 2009. Venoms, venomics, antivenomics. FEBS Lett. 583, 1736e1743. Cardoso, K.C., Da Silva, M.J., et al., 2010. A transcriptomic analysis of gene expression in the venom gland of the snake Bothrops alternatus (urutu). BMC Genomics 11, 605. Casewell, N.R., Harrison, R.A., et al., 2009. Comparative venom gland transcriptome surveys of the saw-scaled vipers (Viperidae: Echis) reveal substantial intra-family gene diversity and novel venom transcripts. BMC Genomics 10, 564. Chang, L., Lin, S., et al., 1999. Genetic organization of alpha-bungarotoxins from Bungarus multicinctus (Taiwan banded krait): evidence showing that the production of alpha-bungarotoxin isotoxins is not derived from edited mRNAs. Nucleic Acids Res. 27, 3970e3975. Chatrath, S.T., Chapeaurouge, A., et al., 2011. Identification of novel proteins from the venom of a cryptic snake Drysdalia coronoides by a combined transcriptomics and proteomics approach. J. Proteome Res. 10, 739e750. Ching, A.T., Paes Leme, A.F., et al., 2012. Venomics profiling of Thamnodynastes strigatus unveils matrix metalloproteinases and other novel proteins recruited to the toxin arsenal of rear-fanged snakes. J. Proteome Res. 11, 1152e1162. Ching, A.T., Rocha, M.M., et al., 2006. Some aspects of the venom proteome of the Colubridae snake Philodryas olfersii revealed from a Duvernoy's (venom) gland transcriptome. FEBS Lett. 580, 4417e4422. Chippaux, J.P., Theakston, R.D., 1987. Epidemiological studies of snake bite in French Guiana. Ann. Trop. Med. Parasitol. 81, 301e304. Chippaux, J.P., Williams, V., et al., 1991. Snake venom variability: methods of study, results and interpretation. Toxicon 29, 1279e1303. Cidade, D.A., Simao, T.A., et al., 2006. Bothrops jararaca venom gland transcriptome: analysis of the gene expression pattern. Toxicon 48, 437e461. Correa-Netto, C., Junqueira-de-Azevedo, I.L., et al., 2011. Snake venomics and venom gland transcriptomic analysis of Brazilian coral snakes, Micrurus altirostris and M. corallinus. J. Proteomics 74, 1795e1809. Creer, S., Chou, W.H., et al., 2002. Offshore insular variation in the diet of the Taiwanese bamboo viper Trimeresurus stejnegeri (Schmidt). Zool. Sci. 19, 907e913. Currier, R.B., Calvete, J.J., et al., 2012. Unusual stability of messenger RNA in snake venom reveals gene expression dynamics of venom replenishment. PLoS One 7, e41888. Daltry, J.C., Wuster, W., et al., 1996. Diet and snake venom evolution. Nature 379, 537e540. De Lucca, F.L., Haddad, A., et al., 1974. Protein synthesis and morphological changes in the secretory epithelium of the venom gland of Crotalus durissus terrificus at different times after manual extraction of venom. Toxicon 12, 361e368. Dolab, J.A., de Roodt, A.R., et al., 2014. Epidemiology of snakebite and use of antivenom in Argentina. Trans. R. Soc. Trop. Med. Hyg. 108, 269e276. Doley, R., Mackessy, S.P., et al., 2009. Role of accelerated segment switch in exons to alter targeting (ASSET) in the molecular evolution of snake venom proteins. BMC Evol. Biol. 9, 146. Doley, R., Pahari, S., et al., 2008a. Accelerated exchange of exon segments in Viperid three-finger toxin genes (Sistrurus catenatus edwardsii; Desert Massasauga). BMC Evol. Biol. 8, 196. Doley, R., Tram, N.N., et al., 2008b. Unusual accelerated rate of deletions and insertions in toxin genes in the venom glands of the pygmy copperhead (Austrelaps labialis) from Kangaroo island. BMC Evol. Biol. 8, 70. Doley, R., Pahari, S., et al., 2010. The gene structure and evolution of ku-wap-fusin (kunitz waprin fusion protein), a novel evolutionary intermediate of the Kunitz serine protease inhibitors and waprins from sistrurus catenatus (Massasauga rattlesnake) venom glands. Open Evol. J 4, 31e41. Droege, M., Hill, B., 2008. The genome sequencer FLX systemelonger reads, more applications, straight forward bioinformatics and more complete data sets. J. Biotechnol. 136, 3e10. Ducancel, F., Guignery-Frelat, G., et al., 1988. Complete amino acid sequence of a PLA2 from the tiger snake Notechis scutatus scutatus as deduced from a complementary DNA. Nucleic Acids Res. 16, 9049. Durban, J., Juarez, P., et al., 2011. Profiling the venom gland transcriptomes of Costa Rican snakes by 454 pyrosequencing. BMC Genomics 12, 259. Durban, J., Perez, A., et al., 2013. Integrated “omics” profiling indicates that miRNAs are modulators of the ontogenetic venom composition shift in the Central American rattlesnake, Crotalus simus simus. BMC Genomics 14, 234.

R.K. Brahma et al. / Toxicon 93 (2015) 1e10 Dutertre, S., Jin, A.H., et al., 2013. Deep venomics reveals the mechanism for expanded peptide diversity in cone snail venom. Mol. Cell. Proteomics 12, 312e329. Eid, J., Fehr, A., et al., 2009. Real-time DNA sequencing from single polymerase molecules. Science 323, 133e138. Fox, J.W., Serrano, S.M., 2005. Structural considerations of the snake venom metalloproteinases, key members of the M12 reprolysin family of metalloproteinases. Toxicon 45, 969e985. Francischetti, I.M., My-Pham, V., et al., 2004. Bitis gabonica (Gaboon viper) snake venom gland: toward a catalog for the full-length transcripts (cDNA) and proteins. Gene 337, 55e69. Fry, B.G., 2005. From genome to “venome”: molecular origin and evolution of the snake venom proteome inferred from phylogenetic analysis of toxin sequences and related body proteins. Genome Res. 15, 403e420. Fry, B.G., Lumsden, N.G., et al., 2003a. Isolation of a neurotoxin (alpha-colubritoxin) from a nonvenomous colubrid: evidence for early origin of venom in snakes. J. Mol. Evol. 57, 446e452. Fry, B.G., Scheib, H., et al., 2012. Novel transcripts in the maxillary venom glands of advanced snakes. Toxicon 59, 696e708. Fry, B.G., Scheib, H., et al., 2008. Evolution of an arsenal: structural and functional diversification of the venom system in the advanced snakes (Caenophidia). Mol. Cell. Proteomics 7, 215e246. Fry, B.G., Wuster, W., 2004. Assembling an arsenal: origin and evolution of the snake venom proteome inferred from phylogenetic analysis of toxin sequences. Mol. Biol. Evol. 21, 870e883. Fry, B.G., Wuster, W., et al., 2003b. Molecular evolution and phylogeny of elapid snake venom three-finger toxins. J. Mol. Evol. 57, 110e129. Gojobori, T., Ikeo, K., 1994. Molecular evolution of serine protease and its inhibitor with special reference to domain evolution. Philos. Trans. R. Soc. Lond B Biol. Sci. 344, 411e415. Gonzalez-Andrade, F., Chippaux, J.P., 2010. Snake bite envenomation in Ecuador. Trans. R. Soc. Trop. Med. Hyg. 104, 588e591. Gubensek, F., Sket, D., et al., 1974. Fractionation of Vipera ammodytes venom and seasonal variation of its composition. Toxicon 12, 167e171. Guo, M., Teng, M., et al., 2005. Crystal structure of the cysteine-rich secretory protein stecrisp reveals that the cysteine-rich domain has a Kþ channel inhibitor-like fold. J. Biol. Chem. 280, 12405e12412. Gutierrez, J.M., Avila, C., et al., 1990. Ontogenetic changes in the venom of the snake Lachesis muta stenophrys (bushmaster) from Costa Rica. Toxicon 28, 419e426. Gutierrez, J.M., Lomonte, B., et al., 2009. Snake venomics and antivenomics: proteomic tools in the design and control of antivenoms for the treatment of snakebite envenoming. J. Proteomics 72, 165e182. Gutierrez, J.M., Theakston, R.D., et al., 2006. Confronting the neglected problem of snake bite envenoming: the need for a global partnership. PLoS Med 3, e150. Harrison, R.A., Ibison, F., et al., 2007. Identification of cDNAs encoding viper venom hyaluronidases: cross-generic sequence conservation of full-length and unusually short variant transcripts. Gene 392, 22e33. Hirabayashi, J., Kusunoki, T., et al., 1991. Complete primary structure of a galactosespecific lectin from the venom of the rattlesnake Crotalus atrox. Homologies with Ca2(þ)-dependent-type lectins. J. Biol. Chem. 266, 2320e2326. Huang, X., Madan, A., 1999. CAP3: a DNA sequence assembly program. Genome Res. 9, 868e877. Jayanthi, G.P., Gowda, T.V., 1988. Geographical variation in India in the composition and lethal potency of Russell's viper (Vipera russelli) venom. Toxicon 26, 257e264. Jia, Y., Cantu, B.A., et al., 2008. Complementary DNA sequencing and identification of mRNAs from the venomous gland of Agkistrodon piscivorus leucostoma. Toxicon 51, 1457e1466. Jiang, Y., Li, Y., et al., 2011. Venom gland transcriptomes of two elapid snakes (Bungarus multicinctus and Naja atra) and evolution of toxin genes. BMC Genomics 12, 1. Jin, A.H., Dutertre, S., et al., 2013. Transcriptomic messiness in the venom duct of Conus miles contributes to conotoxin diversity. Mol. Cell. Proteomics 12, 3824e3833. Johnson, E.K., Kardong, K.V., et al., 1987. Observations on white and yellow venoms from an individual southern Pacific rattlesnake (Crotalus viridis helleri). Toxicon 25, 1169e1180. Junqueira-de-Azevedo, I.L., Ching, A.T., et al., 2006. Lachesis muta (Viperidae) cDNAs reveal diverging pit viper molecules and scaffolds typical of cobra (Elapidae) venoms: implications for snake toxin repertoire evolution. Genetics 173, 877e889. Junqueira-de-Azevedo, I.L., Ho, P.L., 2002. A survey of gene expression and diversity in the venom glands of the pitviper snake Bothrops insularis through the generation of expressed sequence tags (ESTs). Gene 299, 279e291. Kashima, S., Roberto, P.G., et al., 2004. Analysis of Bothrops jararacussu venomous gland transcriptome focusing on structural and functional aspects: Iegene expression profile of highly expressed phospholipases A2. Biochimie 86, 211e219. Kerchove, C.M., Carneiro, S.M., et al., 2004. Stimulation of the alpha-adrenoceptor triggers the venom production cycle in the venom gland of Bothrops jararaca. J. Exp. Biol. 207, 411e416. Kini, R.M., Chan, Y.M., 1999. Accelerated evolution and molecular surface of venom phospholipase A2 enzymes. J. Mol. Evol. 48, 125e132.

9

Kordis, D., Gubensek, F., 2000. Adaptive evolution of animal toxin multigene families. Gene 261, 43e52. Lalloo, D.G., Trevett, A.J., et al., 1996. Neurotoxicity, anticoagulant activity and evidence of rhabdomyolysis in patients bitten by death adders (Acanthophis sp.) in southern Papua New Guinea. QJM 89, 25e35. Lavergne, V., Dutertre, S., et al., 2013. Systematic interrogation of the Conus marmoreus venom duct transcriptome with ConoSorter reveals 158 novel conotoxins and 13 new gene superfamilies. BMC Genomics 14, 708. Leao, L.I., Ho, P.L., et al., 2009. Transcriptomic basis for an antiserum against Micrurus corallinus (coral snake) venom. BMC Genomics 10, 112. Li, M., Fry, B.G., et al., 2005a. Eggs-only diet: its implications for the toxin profile changes and ecology of the marbled sea snake (Aipysurus eydouxii). J. Mol. Evol. 60, 81e89. Li, M., Fry, B.G., et al., 2005b. Putting the brakes on snake venom evolution: the unique molecular evolutionary patterns of Aipysurus eydouxii (Marbled sea snake) phospholipase A2 toxins. Mol. Biol. Evol. 22, 934e941. MacDermot, J., Westgaard, R.H., et al., 1978. beta-Bungarotoxin. Separation of two discrete proteins with different synaptic actions. Biochem. J. 175, 271e279. Mardis, E.R., 2008. Next-generation DNA sequencing methods. Annu. Rev. Genomics Hum. Genet. 9, 387e402. Margres, M.J., Aronow, K., et al., 2013. The venom-gland transcriptome of the eastern coral snake (Micrurus fulvius) reveals high venom complexity in the intragenomic evolution of venoms. BMC Genomics 14, 531. McCue, M.D., 2006. Cost of producing venom in three North American Pitviper species. Copeia 2006, 818e825. Menez, A., 1998. Functional architectures of animal toxins: a clue to drug design? Toxicon 36, 1557e1572. Menezes, M.C., Furtado, M.F., et al., 2006. Sex-based individual variation of snake venom proteome among eighteen Bothrops jararaca siblings. Toxicon 47, 304e312. Minton, S.A., Weinstein, S.A., 1986. Geographic and ontogenic variation in venom of the western diamondback rattlesnake (Crotalus atrox). Toxicon 24, 71e80. Moura-da-Silva, A.M., Paine, M.J., et al., 1995. The molecular cloning of a phospholipase A2 from Bothrops jararacussu snake venom: evolution of venom group II phospholipase A2's may imply gene duplications. J. Mol. Evol. 41, 174e179. Nakashima, K., Ogawa, T., et al., 1993. Accelerated evolution of Trimeresurus flavoviridis venom gland phospholipase A2 isozymes. Proc. Natl. Acad. Sci. U. S. A. 90, 5964e5968. Neiva, M., Arraes, F.B., et al., 2009. Transcriptome analysis of the Amazonian viper Bothrops atrox venom gland using expressed sequence tags (ESTs). Toxicon 53, 427e436. Ogawa, T., Nakashima, K., et al., 1996. Accelerated evolution of snake venom phospholipase A2 isozymes for acquisition of diverse physiological functions. Toxicon 34, 1229e1236. Ohno, M., Menez, R., et al., 1998. Molecular evolution of snake toxins: is the functional diversity of snake toxins associated with a mechanism of accelerated evolution? Prog. Nucleic Acid Res. Mol. Biol. 59, 307e364. OmPraba, G., Chapeaurouge, A., et al., 2010. Identification of a novel family of snake venom proteins Veficolins from Cerberus rynchops using a venom gland transcriptomics and proteomics approach. J. Proteome Res. 9, 1882e1893. Ozsolak, F., Platt, A.R., et al., 2009. Direct RNA sequencing. Nature 461, 814e818. Pahari, S., Bickford, D., et al., 2007a. Expression pattern of three-finger toxin and phospholipase A2 genes in the venom glands of two sea snakes, Lapemis curtus and Acalyptophis peronii: comparison of evolution of these toxins in land snakes, sea kraits and sea snakes. BMC Evol. Biol. 7, 175. Pahari, S., Mackessy, S.P., et al., 2007b. The venom gland transcriptome of the Desert Massasauga rattlesnake (Sistrurus catenatus edwardsii): towards an understanding of venom composition among advanced snakes (Superfamily Colubroidea). BMC Mol. Biol. 8, 115. Pawelek, P.D., Cheah, J., et al., 2000. The structure of L-amino acid oxidase reveals the substrate trajectory into an enantiomerically conserved active site. EMBO J. 19, 4204e4215. Pawlak, J., Mackessy, S.P., et al., 2009. Irditoxin, a novel covalently linked heterodimeric three-finger toxin with high taxon-specific neurotoxicity. FASEB J. 23, 534e545. Pierini, S.V., Warrell, D.A., et al., 1996. High incidence of bites and stings by snakes and other animals among rubber tappers and Amazonian Indians of the Jurua Valley, Acre State, Brazil. Toxicon 34, 225e236. Pla, D., Gutierrez, J.M., et al., 2012. Second generation snake antivenomics: comparing immunoaffinity and immunodepletion protocols. Toxicon 60, 688e699. Pungercar, J., Kordis, D., et al., 1991. Cloning and nucleotide sequence of a cDNA encoding ammodytoxin A, the most toxic phospholipase A2 from the venom of long-nosed viper (Vipera ammodytes). Toxicon 29, 269e273. Pyron, R.A., Burbrink, F.T., et al., 2011. The phylogeny of advanced snakes (Colubroidea), with discovery of a new subfamily and comparison of support methods for likelihood trees. Mol. Phylogenet. Evol. 58, 329e342. Pyron, R.A., Burbrink, F.T., et al., 2013. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evol. Biol. 13, 93. Rajagopalan, N., Pung, Y.F., et al., 2007. Beta-cardiotoxin: a new three-finger toxin from Ophiophagus hannah (king cobra) venom with beta-blocker activity. FASEB J. 21, 3685e3695.

10

R.K. Brahma et al. / Toxicon 93 (2015) 1e10

Ramos-Cerrillo, B., de Roodt, A.R., et al., 2008. Characterization of a new polyvalent antivenom (Antivipmyn Africa) against African vipers and elapids. Toxicon 52, 881e888. Rao, V.S., Swarup, S., et al., 2004. The catalytic subunit of pseutarin C, a group C prothrombin activator from the venom of Pseudonaja textilis, is structurally similar to mammalian blood coagulation factor Xa. Thromb. Haemost. 92, 509e521. Reza, M.A., Minh, L.T., et al., 2006. Molecular evolution caught in action: gene duplication and evolution of molecular isoforms of prothrombin activators in Pseudonaja textilis (brown snake). J. Thromb. Haemost. 4, 1346e1353. Rodrigues, R.S., Boldrini-Franca, J., et al., 2012. Combined snake venomics and venom gland transcriptomic analysis of Bothropoides pauloensis. J. Proteomics 75, 2707e2720. Rokyta, D.R., Lemmon, A.R., et al., 2012. The venom-gland transcriptome of the eastern diamondback rattlesnake (Crotalus adamanteus). BMC Genomics 13, 312. Rokyta, D.R., Wray, K.P., et al., 2011. A high-throughput venom-gland transcriptome for the Eastern Diamondback Rattlesnake (Crotalus adamanteus) and evidence for pervasive positive selection across toxin classes. Toxicon 57, 657e671. Rokyta, D.R., Wray, K.P., et al., 2013. The genesis of an exceptionally lethal venom in the timber rattlesnake (Crotalus horridus) revealed through comparative venom-gland transcriptomics. BMC Genomics 14, 394. Sanger, F., Nicklen, S., et al., 1977. DNA sequencing with chain-terminating inhibitors. Proc. Natl. Acad. Sci. 74, 5463e5467. Schwartz, T.S., Tae, H., et al., 2010. A garter snake transcriptome: pyrosequencing, de novo assembly, and sex-specific differences. BMC Genomics 11, 694. Scott, D.L., Sigler, P.B., 1994. Structure and catalytic mechanism of secretory phospholipases A2. Adv. Protein Chem. 45, 53e88. Serrano, S.M., Maroun, R.C., 2005. Snake venom serine proteinases: sequence homology vs. substrate specificity, a paradox to be solved. Toxicon 45, 1115e1132. Shine, R., 1987. Ecological ramifications of prey size: food habits and reproductive biology of Australian copperhead snakes (Austrelaps, Elapidae). J. Herpetol. 21, 21e28. Siang, A.S., Doley, R., et al., 2010. Transcriptomic analysis of the venom gland of the red-headed krait (Bungarus flaviceps) using expressed sequence tags. BMC Mol. Biol. 11, 24. St Pierre, L., Birrell, G.W., et al., 2007. Diversity of toxic components from the venom of the evolutionarily distinct black whip snake, Demansia vestigiata. J. Proteome Res. 6, 3093e3107. Sunagar, K., Jackson, T.N., et al., 2013. Three-fingered RAVERs: rapid accumulation of variations in exposed residues of snake venom toxins. Toxins. Basel 5, 2172e2208. Tamiya, T., Lamouroux, A., et al., 1985. Cloning and sequence analysis of the cDNA encoding a snake neurotoxin precursor. Biochimie 67, 185e189.

Tsetlin, V., 1999. Snake venom alpha-neurotoxins and other ‘three-finger’ proteins. Eur. J. Biochem. 264, 281e286. Tu, A.T., 1973. Neurotoxins of animal venoms: snakes. Annu. Rev. Biochem. 42, 235e258. Valente, R.H., Guimaraes, P.R., et al., 2009. Bothrops insularis venomics: a proteomic analysis supported by transcriptomic-generated sequence data. J. Proteomics 72, 241e255. Vandenplas, M.L., Vandenplas, S., et al., 1985. Characterization of the messenger RNA population coding for components of viperid snake venom. Toxicon 23, 289e305. Vidal, N., Delmas, A.S., et al., 2007. The phylogeny and classification of caenophidian snakes inferred from seven nuclear protein-coding genes. C. R. Biol. 330, 182e187. Vonk, F.J., Jackson, K., et al., 2011. Snake venom: from fieldwork to the clinic: recent insights into snake biology, together with new technology allowing highthroughput screening of venom, bring new hope for drug discovery. Bioessays 33, 269e279. Voris, H.K., Voris, H.H., 1983. Feeding strategies in marine snakes: an analysis of evolutionary, morphological, behavioral and ecological relationships. Am. Zool. 23, 411e425. Wagstaff, S.C., Harrison, R.A., 2006. Venom gland EST analysis of the saw-scaled viper, Echis ocellatus, reveals novel alpha9beta1 integrin-binding motifs in venom metalloproteinases and a new group of putative toxins, renin-like aspartic proteases. Gene 377, 21e32. Wagstaff, S.C., Sanz, L., et al., 2009. Combined snake venomics and venom gland transcriptomic analysis of the ocellated carpet viper, Echis ocellatus. J. Proteomics 71, 609e623. Wang, Z., Gerstein, M., et al., 2009. RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57e63. Willemse, G.T., Hattingh, J., 1979. Physiological effects of fresh freeze-dried and commercially prepared rinkals (Hemachatus haemachatus) venom. Toxicon 17, 89e93. Willemse, G.T., Hattingh, J., et al., 1979. Changes in composition and protein concentration of puff adder (Bitis arietans) venom due to frequent milking. Toxicon 17, 37e42. Zelanis, A., Andrade-Silva, D., et al., 2012. A transcriptomic view of the proteome variability of newborn and adult Bothrops jararaca snake venoms. PLoS Negl. Trop. Dis. 6, e1554. Zhang, B., Liu, Q., et al., 2006. Transcriptome analysis of Deinagkistrodon acutus venomous gland focusing on cellular structure and functional aspects using expressed sequence tags. BMC Genomics 7, 152.